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We develop a general formalism to describe the propagation of a near-resonant electromagnetic 
field in a medium composed of magnetodielectric resonators. As the size and the spatial separa- 
tion of nanofabricated resonators in a metamaterial array is frequently less than the wavelength, 
we describe them as discrete scatterers, supporting a single mode of current oscillation represented 
by a single dynamic variable. We derive a Lagrangian and Hamiltonian formalism for the coupled 
electromagnetic fields and oscillating currents in the length gauge, obtained by the Power-Zienau- 
Woolley transformation. The response of each resonator to electromagnetic field is then described 
by polarization and magnetization densities that, to the lowest order in a multipole expansion, gen- 
erate electric and magnetic dipole excitations. We derive a closed set of equations for the coherently 
scattered field and normal mode amplitudes of current oscillations of each resonator both within 
the rotating wave approximation, in which case the radiative decay rate is much smaller than the 
resonance frequency, and without such an assumption. The set of equations includes the radiative 
couplings between a discrete set of resonators mediated by the electromagnetic field, fully incorpo- 
rating recurrent scattering processes to all orders. By considering an example of a two-dimensional 
split ring resonator metamaterial array, we show that the system responds cooperatively to near- 
resonant field, exhibiting collective eigenmodes, resonance frequencies, and radiative linewidths that 
result from strong radiative interactions between closely-spaced resonators. 
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I. INTRODUCTION 

Recent advances in nanofabrication provide a variety of 
tools for engineering the electromagnetic (EM) response 
of metamaterials in the radiofrcquency, microwave, and 
optical domains. Metamaterials consist of arrays of ar- 
tificially constructed magnetodielectric resonators which 
typically interact strongly with the incident and scat- 
tered EM fields. These resonator structures frequently 
extend over length scales smaller than the wavelength 
of the EM field with which they interact. For exam- 
ple, a metamaterial might comprise isolated circuit el- 
ements, or meta-atoms, embedded in a dielectric host 
medium. Whereas the EM properties of natural atoms 
are fixed, modifying the design of artificially constructed 
meta-atoms can endow them with a wide range of elec- 
tric and/or magnetic responses. Such control allows one 
to produce materials with EM properties such as nega- 
tive index of refraction^— or negative group velocities' 
These materials could conceivably be employed to create 
perfect lenses^' and electromagnetic cloaks' - — 

The exciting EM phenomena of nanofabricated meta- 
materials can often depend on the effective bulk 
properties of the sample. Homogcnization theories 
have met with substantial success in describing these 
properties'^ - — Homogenization leads to effective con- 
tinuum models that strive to treat excitations using av- 
eraged polarization and magnetization densities formed 
by current oscillations within the unit-cell resonators. 
Analyzing an EM response using uniform medium de- 
scriptions, however, can be complicated by the fact that 



recurrent scattering events, in which a photon scatters 
more than once off the same resonator, produce in- 
teractions which can strongly influence a system's EM 
response^ - — In certain circumstances, the bulk permit- 
tivity and permeability can be inferred by analyzing the 
transmission and reflection properties of a metamaterial 
with finite thickness , 10 ' 11 or from the scattering proper- 
ties of a metamaterial's constituent slabsJ^— But, accu- 
rately accounting for strong interactions between a meta- 
material's unit cells often requires simplifying assump- 
tions such as the elements being arranged in an infinite 
lattice'tr— The discrete translational symmetry of the 
infinite lattice can be exploited, e.g., to approximate the 
local field corrections in a medium of discrete magneto- 
electric scatterers i2i 

The discrete nature of metamaterials becomes appar- 
ent when the infinite lattice symmetry is broken. The 
strongly interacting nature of these structures renders 
them very sensitive to finite size effect o 32 : 33 and to disor- 
der in the latticei 34 ' 35 In systems of discrete resonators, 
interference of different scattering paths between the el- 
ements can result, e.g., in light localization! 36 i 37 This ef- 
fect is analogous to Anderson localization of electrons 
in solids. Even in regular arrays, strong interactions 
between resonators can find important applications in 
metamaterial systems, providing precise control and ma- 
nipulation of EM fields on a subwavelength scale, e.g., by 
localizing sub-diffraction field hot-spots i 38 ' 39 As another 
example, a system of interacting resonant wires was used 
to produce a meta-lens able to transfer subwavelength 
features of an evanescent field to propagating waves 
In essence, recurrent scattering events produce strong in- 
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teractions between meta-atoms that contribute to these 
effects. As a result, ensembles of interacting resonators 
exhibit collective mode of oscillation with discrete reso- 
nance frequencies and radiative emission rates. In princi- 
ple, one can calculate the scattered held profile in a meta- 
material by having knowledge of how the material com- 
prising the meta-atoms reacts to the EM held. One could 
then numerically integrate Maxwell's equations with a 
numerical mesh small enough to resolve the features of 
each meta-atom. This, however, becomes computation- 
ally intractable when the system contains more than a 
few resonator elements. 

In this article, we develop a simplified, computation- 
ally efficient formalism that captures the fundamental 
physical properties of a finite metamatcrial. In partic- 
ular, we show how EM mediated interactions can form 
a cooperative response of the metamaterial's constituent 
resonators. In this model, each unit cell element, or meta- 
molecule, of the metamaterial array is formed by com- 
binations of circuit elements acting as resonators that 
interact with the incident and scattered EM fields. In 
several metamaterial realizations, a metamolccule may 
further be divided into separate sub-elements, e.g., iso- 
lated circuit elements that can naturally be considered 
as the elementary building blocks of the metamaterial 
sample. We refer to such elementary building blocks as 
meta-atoms. We assume each meta-atom supports a sin- 
gle mode of current oscillation represented by a single 
dynamic variable. 

The theoretical formalism we introduce describes the 
collective response of a metamaterial array to an incident 
EM field. To develop this formalism, we begin with the 
Lagrangian and Hamiltonian representations for charge 
and current distributions interacting with EM fields. An- 
alyzing the system in the length gauge, obtained by 
the Power- Zienau-Woolley transformation,— ~— we de- 
rive coupled equations for the EM fields and resonators. 
A single resonator interacting with its self-generated 
fields behaves as an LC circuit in which emission of EM 
radiation damps the current oscillation. An incident EM 
field drives each resonator. However, each meta-atom 
is also driven by fields scattered from all other meta- 
atoms in the metamaterial array. By integrating out the 
EM fields, we derive a set of equations for the meta- 
atom current oscillations which describes the collective 
response of the array to the incident field. Each eigen- 
mode of this system of equations represents a collective 
oscillation distributed over the resonators with a partic- 
ular resonance frequency and radiative decay rate. Some 
modes are superradiant, with emission rates enhanced 
by collective interactions. In other modes, EM mediated 
interactions result in subradiant emission in which ra- 
diation repeatedly scattered between resonators remains 
trapped, slowly leaking away from the metamaterial. As 
an example, we analyze a 2D array of split ring resonators 
and examine several of its characteristic collective modes. 
We find that EM mediated interactions can produce a 
broad distribution of collective emission rates, and that 



the width of this distribution is sensitive to the inter- 
resonator spacing. For example, in a 33 x 33 array in 
which the resonators are separated by half a wavelength 
of the resonant light, the radiative emission rate can be 
suppressed by five orders of magnitude. On the other 
hand, when the spacing is increased to 1.4 wavelengths, 
the emission rate is only suppressed by a factor of five. 

In previous Lagrangian treatments, the interaction be- 
tween elements of a single metamolecule was accounted 
for by a phenomenological coupling between meta-atom 
dynamic variables^ Similar phenomenological coupling 
between nearest neighbor resonators can also describe the 
propagation dynamics of excitations in a one-dimensional 
chain of metamolecules4^ Radiative losses were ac- 
counted for by additional dissipativc terms. However, 
important effects such as superradiance or subradiance 
of collective modes cannot be modeled in this way. By 
contrast, in our treatment the interactions between meta- 
atoms are mediated entirely by the scattered EM fields; 
the radiation lost through decay of one meta-atom can 
drive another and vice-versa. The resulting collective 
modes of the system can, therefore, exhibit cither subra- 
diant or superradiant decay. 

This article is organized as follows. We highlight the 
main results of the developed formalism for the collec- 
tive response of a metamaterial sample to EM fields in 
Sec. |TTJ In Sec. IIII1 wc set up our description of the 
metamaterial. We provide a theoretical description of 
the system dynamics in Sec. IIVI where we introduce the 
Lagrangian and derive the Hamiltonian for the system 
and the equations of motion for the meta-atoms. We 
also arrive at expressions for the scattered EM fields that 
drive the meta-atom dynamics. A derivation of our La- 
grangian from that describing arbitrary charged particles 
interacting with the EM field in the Coulomb gauge is 
provided in Appendix and we elaborate on the deriva- 
tion of the Hamiltonian in Appendix |Bl Wc combine the 
field and meta-atom dynamics in Sec. [V] to arrive at cou- 
pled equations of motion between meta-atoms in the ro- 
tating wave approximation, in which the meta-atom de- 
cay rates are much less than their resonance frequencies. 
A more general model for collective interactions is pro- 
vided in Appendix [C] In Sec. I VII wc apply the theo- 
retical formalism to describe collective modes of oscilla- 
tion in an array of symmetric split ring resonators. Col- 
lective modes of these resonators are connected to the 
linewidth narrowing^ of a transmission resonance ob- 
served in Ref.HU- In Sec. I VIII we quantize this formalism 
in the special case that the resonators do not suffer from 
thermal or ohmic losses. Conclusions follow in Sec. IVIIII 



II. KEY RESULTS: COLLECTIVE DYNAMICS 
ARISING FROM RECURRENT EM 
SCATTERING 

In this section, we summarize key results presented 
in this article. Ultimately, we describe the collective dy- 
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namics arising from an ensemble of magnctodielectric res- 
onators interacting via a near-resonant EM field. When 
such resonators are placed close to each other, the system 
can respond to EM fields cooperatively. In order to pro- 
vide a computationally efficient description, we consider 
a metamaterial array composed of a set of TV discrete 
meta-atoms. We assume each meta-atom j (J = 1 . . . N) 
supports a single mode of current oscillation whose be- 
havior is described by a single dynamic variable Qj , with 
units of charge, and its rate of change Ij = Qj, with units 



of current. As described in Sec. IIII1 the current oscilla- 
tion produces a polarization density Pj (r, t) proportional 
to Qj and a magnetization density Mj(r, f) proportional 
to Ij. An incident wave with electric field E; n (r,i) and 
magnetic induction field Bj n (r,t) impinges on the sys- 
tem. 

For the coupled set of circuit elements and EM 
fields we derive a Lagrangian and Hamiltonian formal- 
ism in Sec. IIV Al The Lagrangian is expressed in the 
length gauge, obtained by the Power-Zicnau-Woollcy 
transformation^^ For the dynamical variables Qj and 
the EM vector potential A(r) we obtain the corre- 
sponding conjugate momenta <j>j [Eq. (JT5J)] and — D(r) 
[Eq. pi]) ], respectively. Here D = eoE + P denotes 
the electric displacement field. The joint dynamics of 
the meta-atom and EM fields are then governed by the 
Hamiltonian 



H = H 



EM 



-/ d 3 r |p| 



E 



— ( d 3 rT>(r,t) •P i (r) 
£o J 



(1) 



where "Hem [Eq. (|32j) ] is the Hamiltonian for the free 
EM field and $j [Eq. p7|) ] is an effective magnetic 
flux through the meta-atom. The final term of the 
Hamiltonian accounts for interactions between electric 
dipoles distributed in the current oscillations and the 
electric field, while magnetic interactions are contained 



$j) and arise in the relationships between 



<&j and Ij. 



From the Hamiltonian we derive a coupled set of equa- 
tions for the EM fields and the meta-atoms. The inci- 
dent EM fields drive current oscillations within the meta- 
atoms, thereby inducing polarization and magnetization 
densities. In Sec. IIV Cl we derive and integrate the equa- 
tions for the total EM fields that are expressed in terms of 
the incident fields and the fields scattered from the polar- 
ization and magnetization densities of the meta-atoms. 
Specifically, currents in meta-atom j, when oscillating 
at a frequency CI, produce the monochromatic scattered 



electric field Egj and magnetic field Hs,j given by 



Es,i(r,fi) 



47re n 



G(r-r',Q)-P 3 -(r',Q) 



f 



G x (r-r',n)-M j (r',Cl) , (2a) 



H S j(r,f2) = — / dV G(r — r', Cl) ■ Mj(r', SI) 

-cGx(r-r , ,f2).P j (r / J fi)] > (2b) 

where G(r — r', CI) is the radiation kernel connecting an 
oscillating electric (magnetic) dipole at position r' to the 
electric (magnetic) field at position r, while G x (r — r', CI) 
connects an electric (magnetic) dipole at r' to its radi- 
ated magnetic (electric) field at r4i Expressions for these 
radiation kernels are given in Eqs. (|64|) and (|65p . 

The total electric and magnetic fields are obtained as 
a sum of the incident fields and the fields scattered by all 
the meta-atoms in the system 

D(r,t) =D in (r,i) + £D SJ (r,i), (3) 

3 

B(r,i) =B in (r,i) + ^Bsj(r,i), (4) 



where we have the scattered magnetic induction Bsj = 
/io(Hs,j + Mj) and the scattered electric displacement 
Dsj = £oEs,j + Pj from the meta-atom j. 

Although, according to Eqs. ([2|), the polarization and 
magnetization densities of all the meta-atoms act as 
source terms that determine the scattered EM fields, 
there is, in general, no simple way of solving for the po- 
larization Pj(r, CI) and magnetization Mj(r, CI) densities 
themselves. The equations for near-resonant EM fields 
and closely-spaced resonators are strongly coupled, and 
the meta-atoms are driven by both the incident fields and 
fields scattered by all other meta-atoms in the system. 
This is illustrated by Hamilton's equations of motion for 
the resonators 



Qj(t) = Ij(t), 



(5a) 
(5b) 



where the total electric field induces an effective elec- 
tromotive force (EMF) £j [Eq. (fTT)) ]. driving the meta- 
atoms. 

Solving the coupled dynamical equations for the res- 
onators ([5|) and the EM fields ([2|) constitute the central 
results of the paper. We begin in Sec. IV Al by consider- 
ing a single meta-atom. A meta-atom not only feels the 
influence of the incident EM field and the fields scattered 
from other meta-atoms, but also its self-generated field. 
We show that interactions between a meta-atom j and 
its self-generated field produces an effective damped LC 
circuit for the current oscillations with capacitance Cj 
[Eq. self-inductance Lj [Eq. ([55)1 ] and resonance 

frequency uij = 1/ J Cj Lj . The oscillating electric and 
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magnetic dipoles of the meta-atom scatter EM fields and 
therefore induce a radiative decay at rates Tej and Tmj, 
respectively. 

In a metamaterial array of several meta-atoms we then 
solve the coupled set of equations ([5]) and ([2]) when each 
meta-atom is also driven by the scattered fields from all 
the other meta-atoms. This results in multiple scattering 
events and yields EM field mediated interactions between 
the meta-atoms. In particular, when the multiple scat- 
tering between closely-spaced resonators becomes domi- 
nant, so that the EM wave is scattered more than once 
by the same scatterer (this is called recurrent scatter- 
ing), the system responds to EM fields cooperatively. In 
order to analyze the eigenmodes of such a system, it is 
beneficial to introduce the excitation amplitudes of each 
meta-atom LC circuit in terms of its dynamical coordi- 
nates and the canonical momenta. In particular, when 
the decay rates are much less than the resonance frequen- 
cies, the meta-atom dynamics are well described by the 
slowly varying normal variables 



b 3 (t) 




(6) 



In terms of the normal variables we then derive a linear 
set of equations for the meta-atoms whose interactions 
are mediated by the EM fields by explicitly integrating 
out the scattered fields 



(7) 



The expressions for the coupling matrix C between the 
meta-atoms and the driving terms by the incident fields 
fin.j are derived in Sec. IV Bl The diagonal elements of C 
reflect the resonance frequencies and decay rates of the 
single meta-atoms in isolation, while the off-diagonal ele- 
ments arise from scattered electric and magnetic fields 
interacting with the meta-atom electric and magnetic 
dipoles. The coupled equations include the recurrent 
scattering events between the meta-atoms to all orders. 
We generalize our treatment to account for stronger in- 
teractions between meta-atoms, i.e., where interactions 
mediated by scattered fields are comparable to the ef- 
fects of the self-generated field, in Appendix [CJ 

A system of N meta-atoms supports N collective 
modes of current oscillation, each matched to an eigen- 
vector of the matrix C. Each mode i has its own collective 
resonance frequency and decay rate given in terms of its 
eigenvalue Aj as 



Qi = -Im(Ai) + fi , 
7i = -2He(Ai), 



(8a) 
(8b) 



respectively. As a result of the interactions, the collec- 
tive emission rates can be either much less than (subradi- 
ant) or much greater than (supcrradiant) the constituent 
single meta-atom decay rates. We demonstrate this in 
Sec. I VII where we consider the collective effects on a 2D 



metamaterial array of symmetric split ring resonators 
(SRRs), metamolecules possessing reflection symmetry 
which consist of two concentric circular arcs of equal 
length. Even in a relative small metamaterial sample 
of 33 x 33 unit-cell resonators for the lattice spacing of a 
half-wavelength we find that the spectrum of resonance 
frequencies exhibits a long tail of strongly subradiant 
eigenmodes. The most subradiant mode of the system 
possesses a radiative decay rate of about five orders of 
magnitude less than that of an isolated meta-atom. This 
eigenmode exhibits a checkerboard phase-pattern of dom- 
inantly electric dipole excitations. We also find that the 
strong response of the metamaterial sample is very sen- 
sitive to the spacing between the resonators. We analyze 
the spectrum for the lattice spacing of 1.4 wavelength in 
which case the distribution of the decay rates is consid- 
erably narrowed. The most subradiant mode now has a 
resonance linewidth that is five times narrower than the 
one of the isolated unit-cell resonator. Finally, we also 
provide an example how the propagation dynamics of ex- 
citations in a metamaterial array can be analyzed using 
the collective eigenmodes. We find that the lattice spac- 
ing, and hence the interactions between the resonators, 
strongly influence the rate at which excitations spread 
over the array. 



III. DISCRETE RESONATOR MODEL OF A 
METAMATERIAL 

To develop the formalism characterizing interactions of 
magnetodielectric resonators in EM fields, we first pro- 
vide a detailed description of the metamaterial and the 
model we use to represent it. We consider an ensem- 
ble of metamolecules, unit-cell elements that comprise 
the metamaterial, driven by an incident EM field. Each 
metamolecule can be decomposed into some number of 
meta-atoms, which may correspond, for example, to in- 
dividual circuit elements. We model our metamaterial 
as an ensemble of N meta-atoms. The position of the 
meta-atom j is denoted by Rj (j = 1, . . . , N). An ex- 
ternal beam with electric field Ej n (r, t) and associated 
magnetic induction Bj n (r, f) impinges on the ensemble, 
driving the meta-atoms. We assume the incident field is 
bandwidth limited with a spectrum centered at angular 
frequency S7o, and that the spatial extent of each meta- 
atom lies well within a carrier wavelength A = 2itc/Qq. 

The meta-atoms may be composed of, e.g., metallic 
circuit elements supporting plasmonic oscillations, allow- 
ing charges and currents to flow internally. The cur- 
rent and charge distributions produce EM fields, which 
in turn, influence the dynamics of these distributions. 
As such, each element supports various eigenmodes of 
current oscillation ! 48 ' 49 For simplicity, we identify meta- 
atom j with a single eigenmode of current oscillation 
whose state can be described by a single dynamic vari- 
able Qj(t) with units of charge and whose spatial pro- 
file is described by time-independent functions Pj (r) and 
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Wj(r). These mode functions are denned such that the 
polarization Pj(r, t) and magnetization M : (r, t) densi- 
ties associated with atom j are 

P j (v,t)=Q j (t)p j (v), (9a) 
M i (r J t)=J J (f)w J -(r), (9b) 

where Ij(t) = dQj/dt is the current. The definitions 
of the polarization and magnetization lead to the ex- 
pressions of the charge and current densities within each 
meta-atom, 

p j {r,t) = -Q j (t)V-p j (r) (10a) 
j^r.t) = ( Pi (r)+Vx Wj (r)) I s (t). (10b) 

The total polarization and magnetization densities of 
the system are obtained from a sum over the polariza- 
tion and magnetization densities of every meta-atom j, 
respectively, 

P(r,t)=^P i (r,t), (11) 

3 

M(r,t) = ^M i (r,t). (12) 

3 

We choose the mode functions so that they are zero out- 
side the neighborhood of the meta-atom. 

Wc note that, in general, the various parts of a circuit 
element contain charge and current densities that could 
behave independently of one and other; they could there- 
fore be represented by separate dynamic variables. These 
extra degrees of freedom could be described by assigning 
multiple modes of current oscillation to the element, each 
with its own dynamic variable and mode functions to de- 
scribe the corresponding polarization and magnetization 
densities. The resulting set of mode function dynamic 
variables could then interact with one and other via the 
EM fields. In essence, one could view a circuit element 
as an ensemble of meta-atoms that touch or overlap with 
one and other. In this work, however, we have assumed 
that the mode functions have been chosen so that they 
are eigenmodes of elements, i.e., there is a zero net in- 
teraction between the modes in a given circuit element 
within a metamolecule. We therefore identify a meta- 
atom with an eigenmode of current oscillation within a 
circuit element and treat each meta-atom as possessing 
only a single mode of current oscillation. This is analo- 
gous to approximating an atom interacting with the EM 
field as a two-level atom. In the present work, we will 
not address how the eigenmodes of current oscillations 
are calculated. For isolated circuit elements they could 
be computed numerically solving Maxwell's equations us- 
ing actual material parameters. Alternatively, one could 
obtain the meta-atom resonance properties directly from 
experimental measurements, or estimate them using ge- 
ometrical arguments. 



IV. SYSTEM DYNAMICS 

In this section we introduce the Lagrangian and Hamil- 
tonian formalism for a magnetodielectric medium inter- 
acting with EM fields, specifically derived for a system 
consisting of circuit elements whose dynamic variables 
represent eigenmodes of current oscillations. The La- 
grangian is expressed in the length gauge, obtained by the 
Power-Zienau-Woolley transformation^^— We find that 
this particular representation of electromagnetism turns 
out to be especially useful for describing localized, collec- 
tively interacting circuit elements. The specific details of 
the Power-Zienau-Woolley transformation are covered in 
Appendix [A] 

From the Lagrangian we derive the conjugate momenta 
for the dynamic variables of the meta-atoms and the EM 
fields, and the Hamiltonian for the system. The dynam- 
ics of the model describe charge and current densities of 
the system interacting with the EM fields. Wc derive a 
coupled set of equations for the EM fields and the res- 
onators in which both the electric and magnetic fields 
drive the meta-atom dynamics. The expressions for the 
electric and magnetic fields are obtained in terms of the 
incident fields illuminating the sample and the fields scat- 
tered from the polarization and magnetization densities 
that represent the meta-atoms in the medium. 



A. The Lagrangian and Hamiltonian formalism for 
meta-atoms interacting with EM fields 

We treat the dynamics of the system in the Coulomb 
gauge beginning with the Lagrangian formalism. It 
is particularly advantageous to study the EM response 
in a gauge representation obtained by the Power- 
Zienau-Woolley transformation^^— In Appendix El we 
show that the Lagrangian in the Power-Zienau-Woolley 
picture^ can be written in terms of meta-atom dynamic 
variables as 

C = K + £em + Vboui + J2 IQi $) £ i + h ' ( 13 ) 

3 

where K, is an effective kinetic energy given by 

K = £^. (14) 

3 

The phenomenological kinetic inductance lj of meta- 
atom j provides, within the effective single- particle de- 
scription of the system, an inertia to the current oscil- 
lation that would be present in the absence of EM field 
interactions. This inertia could result, for example, from 
the effective mass of charge carriers or surface plasmons 
within the meta-atom. Excitation of a current oscillation 
displaces charge carriers from their equilibrium configu- 
ration producing a charge density Pj(v) [Eq. (|10ap ] within 
meta-atom j. The meta-atom charge densities interact 
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via the instantaneous Coulomb interaction 

*i >s / P;( r )ft>'( r 



(15) 



The current oscillation in meta-atom j interacts with 
external EM fields via an effective electromotive force 
(EMF) £j and an effective magnetic flux $j(£) through 
that meta-atom: 



£j(t) = J d 3 rE(r,i). Pj (r), 
^j(t) = J d 3 rB{r,t) ■ Wj(r) . 



(16) 
(17) 



The EMF interacts with the charge Qj on the circuit, 
while the current Ij interacts with the magnetic flux. 
The Lagrangian for the free EM field £em is given in 
terms of the Coulomb gauge vector potential A as 



C 



EM 



£0 
2 



<Tr 



dA 
~dt 



c 2 IV x A| 2 



(18) 



The Lagrangian for the free EM field represents the ra- 
diative fields that are responsible for the excitations of 
the meta-atoms. 

We now wish to determine the Hamiltonian for the sys- 
tem. We proceed by identifying the conjugate momenta 
of the dynamic variables. Those for charges are given by 



dl~ 



(19) 



Note that in the limit where lj is vanishingly small, the 
conjugate momentum of the charge is dominated by the 
flux through the circuit. This is often the case in mi- 
crowave metamaterials, where the EM interactions dwarf 
the effects of charge carrier inertia. The vector poten- 
tial represents a continuous field of dynamic variables 
which possess a corresponding continuum of conjugate 
momenta defined as 



n(r,t) = 



dC 



dA(r,t) 



(20) 



This conjugate momentum will have a contribution from 
Cem and pick up a contribution from the interaction term 
J2j£jQj = J d 3 rP ■ E. For a system of neutral meta- 
atoms, the conjugate momentum of the vector potential 
is given by££ 



where 



n(r,t) = -D(r,t) ; 



D = e E + P 



(21) 



(22) 



is the electric displacement field. 

In treating the field dynamics, it is often convenient 
to express these fields in terms of the normal variables 



flq,A(£) describing a plane wave with wavevector q and 
transverse polarization e q A. These normal variables are 
defined such that the electric displacement and magnetic 
fields are given by 



D(r,t)=« f rf 3 (? e g ^eq,Aa q ,A(t)e lq r + C. 
B(r,t)=iM f d 3 ^ g ]T(qxe q , A ) 



(23) 



aq.\{t)e 



iq-r 



+ C.c, 
respectively, where 



2(2tt) ; 



(24) 



(25) 



The normal variables for the EM field satisfy the follow- 
ing relations in terms of the Poisson brackets^ 

{a q ,A, a q\A'} = -iSx,X'S(q- q'), (26) 

and {a q! A, Oq'.A'} = {«q,A, Qj} = {a q .A, 4>j} = 0. 

Having obtained the conjugate momenta and normal 
field variables, one may write the Hamiltonian for the 
system by applying the Legendre transform 



U = STQ j <t>j + j 



d 3 rk-U-C. 



(27) 



It is beneficial to decompose the Hamiltonian H = 
H mra + He into a component containing contributions 
from the meta-atom conjugate momenta, W mm , and a 
component accounting for electric field interactions and 
the free EM field, He- The former contribution is given 
explicitly by 



H, 



E 1 



K. 



But, because <f>j - $j = II j [Eq. ([19])]. 



H, 



21, 



(28) 



(29) 



reduces to the kinetic energy of the current oscillations. 
The terms involving the electric field contribution, on the 
other hand, are given explicitly by 



d r A • D + E • P 



Vco 



C 



EM 



(30) 



It is beneficial to simplify the contribution of Eq. (|30| . 
We carry out this simplification in Appendix [Bj The to- 
tal system Hamiltonian may thus be written in the the 
Power-Zicnau-Woolley picture as^ 



H = Hem 



1 

2T 



d:''r\V\ 



E 



i 

2L 



— / d 3 rT>(r,t) ■ Pj(r) 
£o J 



(31) 
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where the Hamiltonian for the free EM field is 



Hem =~ / 



IV x A| 



(32) 



To understand the dynamics that will arise from this 
Hamiltonian, we examine the physical role of each term 
individually. The interaction between the displacement 
field and polarization density can be written in terms of 
the emitter dynamic variables as 

. f . PjM) = _9i I D(r,t) • Pj -(r) . 

(33) 

This represents an interaction energy between the elec- 
tric displacement and the spatial distribution of electric 
dipoles contained in the polarization density. On the 
other hand, the interaction with the magnetic field be- 
comes apparent when expanding K., which yields 



21, 



21 ' 07 _ ' 



21, 



(34) 



The interaction of meta-atom j with the magnetic field 
arises in the second term. The physical significance 
of this interaction can be understood by expressing 
that contribution in terms of the magnetization density 
[Eq. ©] as 



(i 3 rM., B ^-1. 



3 

T, 



(35) 



Equation (|35|) effectively contains the interaction be- 
tween the magnetization density and the magnetic field. 
Additionally, Eq. (|35|) includes a term proportional to 
the square of the magnetic flux. This artifact appears 
because the magnetization density is a function of the 
meta-atom current Ij rather than its conjugate momen- 
tum. When this portion of the interaction is written as 
entirely in terms of the meta-atom conjugate momen- 
tum, the term proportional to the square of the flux 
disappears. The last term of Eq. (|34p represents a dia- 
magnetic interaction proportional to the square of the 
magnetic field flux through a meta-atom. These inter- 
actions are analogous to the effective magnetization and 
diamagnetic interactions found in the Hamiltonian for 
electrically charged point particles in Ref. [50L 

Finally, we examine the local polarization self-energy 
term appearing in the Hamiltonian [Eq. (|3"2"|) ]. This can 
be expressed in terms of the dynamic variables as 



^-J rf 3 rP(M)-P(r,i) 

= ^^| d 3 rp . (r) . p ., (r) 



(36) 



When the meta-atoms are spatially separated, however, 
their polarization mode functions do not overlap, i.e., 
Pj ■ pji = for j 7^ j'. The presence of P(r, t) ■ P(r, t) 
results only in an interaction of the meta-atom with itself, 
which manifests itself as 

±-jd 3 r-p(r,t)--p(T,t) = J2^- Q fd 3 r | Pj (r)| 2 . 

(37) 

If the meta-atoms were to overlap, a contact potential 
proportional to QjQj> would appear between the over- 
lapping elements j and j' . In the initial Lagrangian 
[Eq. (|13p]. direct interactions between the meta-atoms 
appeared via the instantaneous Coulomb interaction. An 
advantage of the Hamiltonian treatment in the Power- 
Zienau-Woolley picture is that such interactions do not 
appear explicitly; other instantaneous, non-causal contri- 
butions to the dynamics cancel out those of the Coulomb 
interaction. This leaves the meta-atom dynamic vari- 
ables to interact exclusively with the vector potential and 
its conjugate momentum. Any effective interactions be- 
tween meta-atoms are thus mediated by these field dy- 
namic variables^ 



B. The meta-atom dynamics 

The meta-atoms' interaction with the EM fields are 
illustrated by Hamilton's equations of motion describing 
the current oscillations 



Q j{t) = {Qj ,H}=i j (t ) = ^m 

fa{t) = {<t> J ,H}=£ j (t). 



(38a) 
(38b) 



The conjugate momentum is driven entirely by the 
EMF £j, while Eq. (|38a[) is nothing more than a state- 
ment that the rate of change of the charge is the current. 
At first glance, it may appear that the magnetic field docs 
not drive the meta-atoms. However, its effects manifest 
themselves indirectly through a relationship between the 
conjugate momentum <fij and the current Ij that will be 
discussed in Sec. [Vj Effective interactions between the 
resonators come about through multiple scattering of the 
EM field between resonators. 



C. The scattered EM fields 

In the previous subsection we derived the equations for 
the meta-atoms driven by EM fields. In order to find a 
coupled set of equations for the fields and the resonators, 
we need find how the EM fields depend on the state of 
the meta-atom charges and currents. In this subsection 
we derive integral expressions for the scattered EM fields 
where the metamaterial medium acts as a source with 
effective polarization and magnetization densities. The 
total electric and magnetic fields are then represented as 



8 



sums of the incident fields and the scattered fields from 
the medium. The resulting equations for the resonators 
and the EM fields are strongly coupled: the resonator dy- 
namics are driven by the EM fields and the fields them- 
selves depend on the excited meta-atom current oscilla- 
tions. 

We begin with the equation of motion for the normal 
field operators 



da, 



dt 



JcqtSl f d 3 r / g-tq-r' 
£0 J 



(e q , A -P(r',i)) 



+ -(qxe q , A )-M(r',t) 

c 



(39) 



component, with |kj n | = CIq/c, 



D in (r,t) 



= i(k in -r-n t) 



C.c. 



B in (r,t) = W— k in x Din(r,i). 



(45) 
(46) 



One obtains explicit expressions for D and B, by sub- 
stituting Eq. (|30"j) into Eqs. (|23[) and IpH)) . summing over 
the two transverse polarizations e q A , and integrating 
over q. Following this procedure, one obtains the scat- 
tered fields 

D s (r,t) = f I dV[S(r-r',t-t')-P(r',0 

J — OO J 

+ is x (r-r',t-t')'M(r',t')] (47) 



where the first term results in the oscillation of the free 
EM field and in the second term we find the polarization 
and magnetization densities arising from meta-atom cur- 
rents that act as sources for radiation. Upon integrating 
Eq. ([55)1 . one obtains 



& -ic q t (in) 
q,A 



d 3 r' / dt' e"^*-*' V iq ' r ' 



6q,A -P(rV 



+ -(qxe q , A )-M(r',t') 

c 



(40) 



where = limt n _ > _ 00 e lcqt °aq.\(to) is the initial state 



*q A — milt _>_oo c. - u, q ., 

of the plane wave normal variable before it interacts with 
the resonators. 

The incident fields (Di n ,Bj n ) and the scattered fields 
(Dg, Bg) radiated by the meta-atoms comprise the total 
electric displacement and magnetic induction fields 

D(r ) t)=D in (r J t)+D s (r,t), (41) 
B(r,t)=B in (r,t)+Bs(r,t), (42) 

D s (r,t)=^D SJ (M) ; (43) 

3 

B s (r,t)=^B SJ (r,t), (44) 



where Ds,j(r, t) and Bsj(r,i) denote the fields emitted 
by the meta-atom j. The incident fields Di n and B; n are 
given in the plane- wave representation by Eqs. (f23f and 

dH]) with e- lcqt a^l substituted for a q , A . For the field ob- 
servation point r located outside the meta-atoms we have 
Dj n (r, t) = eoEi n (r, t) from Eq. (|22|) . A common situa- 
tion in experiments corresponds to an illumination of a 
metamatcrial sample by a non-focused, monochromatic 
beam that can be approximated by a single plane-wave 



and 
B s (r,t) 



Mo 



[ d 3 r' [S(r-r',t-t')-M(r',t') 

3G J 

-cS x (r-r',t-t')■P(r' 7 t , )], (48) 



where S is the propagator that takes the electric (mag- 
netic) field from the electric (magnetic) dipole source r' 
to the observation point at r, and S x represents the 
propagation of the radiated electric (magnetic) field from 
the magnetic (electric) dipole sources to the observation 
points. The propagator S is given by 



S(r,t) 



d 3 qfc(l-qq) e iq r (e 



-icqt 



Jcqt\ 



IC 

16n~ 3 _ 

^(yy^^) 5{r - ct) - 5{r + d) . (49) 
4-7T r 



The two delta-functions produce retarded and advanced 
time contributions to the scattered fields, with the ad- 
vanced time's contribution arising only at r = 0. The 
derivatives acting on 1/r result in a contact interaction 
proportional to S(r). At first glance, the retarded and ad- 
vanced time contributions may appear to cancel out at 
r = 0, thus nullifying such a contact interaction. How- 
ever, by examining the frequency components of S, one 
can show that this contact interaction does survived 
The corresponding expressions for S x arc 

S x (r,<) = J dV'^^' + e^qxl 

= -V x -!— % [5(r - ct) - 6(r + ct)} (50) 
47rr at 



From the oscillator equations of motion [Eqs. (1381) ]. 
we find that the fields emitted from one resonator will 
drive all of the others. The driven resonators, in turn, 
re-scatter these fields to yet other resonators in the meta- 
material. To more easily account for the cumulative ef- 
fects of these multiple scatterings and identify collective 
modes in the system, we analyze the field and oscillator 
dynamics in the frequency domain. 
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Wc therefore decompose the source fields P and M into 
their frequency components and compute the scattered 
monochromatic constituents of the EM fields. Specifi- 
cally, we write an arbitrary source field, 



f(r,t) 



1 



dOf(r,0)e 



-iflt 



(51) 



in terms of the Fourier components f (r, O). In evaluating 
the response of the field to each monochromatic source 
component, one encounters integrals of the form 



J — C 



dt' S(r - r', t - t')f (r', Sl) e - int ' 

= e -' lf2t S(r-r',0)f(r',0) (52) 



and 



( dt'S x (r-r',t-t')f(r',n)e- 

J — oo 



iflt' 



-iUt 



S x (r-r',0)f(r',0), (53) 



where S(r,0) and S x (r, O) are the monochromatic ver- 
sions of the expressions given in Eqs. (|4"9")l and ([50]) that 
describe the propagation of the radiated fields from the 
source to an observation point. In evaluating these prop- 
agators, we find it convenient to treat positive and neg- 
ative frequencies O separately. We therefore decompose 
the propagators as 

S(r,0) = S (+) (r,ft)8(0) + S (_) (r, 0)6(-0) (54) 
S x (r,0) = S ( x +) (r,ft)6(ft) + S ( x _) (r, fi)6(-0) (55) 

where is the Heaviside function, and the propagators' 
positive and negative frequency components are given by 



and, 



S (±) (r,0) = — (VV-1V 2 ) — 

47r r 



(56) 



(57) 



where k = |0|/c is the angular wavenumber of the radia- 
tion emitted from a monochromatic source of frequency 

o. 

One of our goals is to provide radiated electric and 
magnetic fields E and H, respectively. These are related 
to the electric displacement D and magnetic induction B 
by the familiar expressions 

E(r,t) = l[D(r,t)-P(r,i)], (58a) 

H(r, t) = — B(r, t) - M(r, t) . (58b) 
Mo 



We thus define dimensionlcss radiation kernels 

47r r- 



G(r,0) 



k 3 
An 



S(r,0)-tf(r) 



G x (r,0) = -^S x (r,fi), 



(59) 
(60) 



where, by the relations of Eqs. (|58ap and ()58b[) . the delta 
function in Eq. (fSU)) transforms the monochromatic prop- 
agators of D and B to those of E and H, respectively. 
The Fourier components of the corresponding EM fields 
are thus given by 



E(r, O) = -D in (r, O) + V E s .,(r, O) 

3 

(r, SI) = — B in (r, Sl) + Y, H Sj (r, SI) 



H 



/'o 



(61) 
(62) 



where the fields scattered from meta-atom j are 
k 3 



E S ,,(r,0) 



d 3 r' 



G(r - r', SI) -Pjir', SI) 



1 



-G x (r-r',0) -Mjir^Sl) , (63a 



H S>J -(r,0) = 



4tt 



j3„/ 



G(T-r',Sl)-Mj(r',Sl) 



-cG x (r-r',Sl)-P j (r',Sl) 



(63b) 



As with the monochromatic propagators, we decompose 
the radiation kernels into their positive and negative fre- 
quency components as 

G(r, SI) = G (+) (r, 0)6(0) + G ( " } (r, 0)6(-0) . (64) 

G(r,0) = G ( x +) (r, 0)9(0) + G ( > 7 ) (r,0)e(-0). (65) 

The radiation kernel G^(r — r', O) corresponds to a fa- 
miliar expression for a radiated electric (magnetic) field 
at the observation point r, originating from an oscillat- 
ing electric (magnetic) dipolc residing at r'.— Similarly, 
an oscillating electric (magnetic) dipole at r' generates a 
magnetic (electric) field at r that is represented by the 
radiation kernel G x (r — r', O). The explicit expressions 
for these read 



G (±) (r,0) 



4tt 



' : lh^\kr) 



I6(kr) . 



rr 1 



l)h?\kr) 



j p ±ikr 

G x ± )(r ) 0) = ± ivxV 



(66) 
(67) 



where hn and hi, ' are the spherical Hankel functions 
with order n of the first and second kinds, respectively, 
defined by 

(68a) 



hf ) {x)=±i 



1 



± i- 



(68b) 
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Equations (j6"Tj) - (f6"5"j) . together with the radiation ker- 
nels of Eqs. (|66p and (j6"?| . constitute the main results 
of this subsection. They provide the total electric and 
magnetic fields both inside and outside the metamate- 
rial sample as a function of polarization and magnetiza- 
tion densities that are produced by oscillating currents 
in the meta-atoms. Although we have derived the inte- 
gral expressions for the scattered EM fields in terms of 
the resonator excitations in Eq. (|63|) . in general there is 
no simple way of solving for P(r) and M(r). Together 
with Hamilton's equations for the dynamic variables of 
the electric charges of the meta-atoms [Eqs. (|38|) ]. the 
formulas for the radiated fields form a coupled set of 
equations for the EM fields and the matter. The scat- 
tered fields from each meta-atom drive the dynamics of 
the other meta-atoms in the system, with the EM fields 
mediating interactions between the resonators. For the 
case of near-resonant field excitation and closely-spaced 
circuit elements the coupling between the EM fields and 
the meta-atoms can be strong due to multiple scattering 
processes leading to collective behavior of the system. 

In evaluating the scattered fields of Eq. (|63|). we note 
that because /4 ( r ) contains a 1/r 3 divergence near 
r = 0, the spatial integral of G^(r, 0) in Eq. (|66|) is 
not absolutely convergent around the origin. However, 
as Ref. [23| points out one can handle such a singularity 
by carving an infinitesimal spherical region around r = 
from the integral and treating this region separately. The 
integral over the radiation kernel (|66j) is then defined us- 
ing the convention that the term inside the brackets van- 
ishes over an infinitesimal volume enclosing the origin. 
Mathematically, this is achieved by carrying out the inte- 
gral in this region in spherical coordinates, first integrat- 
ing over the spherical angles, so that only the ^-function 
contributes to the integral.— With this integration pro- 
cedure, the 5-function appearing in Eq. (|66|) is required 
for the scattered fields to satisfy Gauss' law, as well as 
to produce the correct Maxwell's equations, V • D = 
for a neutral system and V • B = 0. The requirement 
that these conditions arc satisfied also confirms that we 
have duly selected the correct field terms in the Hamilto- 
nian (|3Tj) (e.g., electric displacement, instead of electric 
field) and that the integration procedure of the contact 
terms [Eq. (|49| ] has been performed correctly. While the 
(5-function singularity in G does not play a role in the in- 
teractions between non-overlapping meta-atoms, we find 
in Sec. IV Al that it does contribute to interactions of a 
meta-atom with its self-generated field. 



The EM fields derived from the Hamiltonian are indeed 
consistent with Maxwell's equations. To verify this, we 
check that the positive and negative frequency compo- 
nents of a monochromatic field with wavenumber k sat- 



isfy the wave equations with sources and M^— : 
(V 2 + jfc 2 )D^ = -Vx(VxP (± >) 

Tt-V x (69) 
c 

(V 2 + k 2 )B^ = -/i V x (V x M^) 

±i> cfcV x P (±) (70) 

We confirm that the total fields produced by our system 
satisfy Eqs. (f6"9")) and (fTUj) by applying the operator (V 2 + 
fc 2 ) to the total electric and magnetic fields [Eqs. (f6Tj) 
and (|62j) ] . Because the incident waves are composed of 
superpositions of plane waves, the action of the operator 
(V 2 + fc 2 ) on these fields trivially reduces to 

(V 2 + fc 2 )D in = (V 2 + fc 2 )B in = 0. (71) 

Therefore, the only contributions to (V 2 + fc 2 )E and 
(V 2 + fc 2 )H come from the scattered fields Esj and Hsj 
[Eqs. (|63a[) and (|63b|) ] . These contributions are most 
readily determined by expressing the tensor components 
of the radiation kernels in the differential form 



G 



i d 
k^d^ 



kr 



k 2 

47r<5 MT ,(5(fcr 



kr 



(72) 



(73) 



Because the differential operators involved in the radi- 
ation kernels readily commute with (V 2 + fc 2 ), the ex- 
pressions for this operator acting on the scattered fields 
involve contributions of the form 



(V 2 



±ifc r-r' 



r - r 



(74) 



appearing under the integral. Physically, the (5-function 
represents a point source away from which a monochro- 
matic spherical wave {e ±lkr /r) propagates. The resulting 

expressions for (V 2 + fc 2 )E^ and (V 2 + fc 2 )H^ thus 
contain integrals over J-functions which are readily eval- 
uated. Explicitly, for the component fi of the scattered 
electric field, we have 



fc 2 ) Z>gf 



vMp (±) ' 



.fc d 



M 



±)r, 



where 



(75) 



(76) 



is the fj, component of the scattered displacement field 
from meta-atom j. Adding the contributions of Eq. ([75)) 
for all meta-atoms j, produces the equivalent of the 
wave equation [Eq. ([69|], which is the desired result. 
Similarly, one finds that by adding the contributions 
^V(V 2 + fc 2 )i?g^ M for all meta-atoms, one recovers the 
wave equation for the magnetic field [Eq. (|70p]. 
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V. META-ATOM INTERACTIONS MEDIATED 
BY THE EM FIELD 

In the previous section we established how current os- 
cillations in the meta-atoms respond to the EM field 
[Eqs. ([38]) ]. Additionally, we arrived at expressions for 
the electric and magnetic fields scattered by the meta- 
atoms [Eqs. (|6Tj) - (|65l) ]. These fields were solved in terms 
of the magnetization and the polarization densities, gen- 
erated by the resonator excitations. The current oscilla- 
tions in each meta-atom thus depend on the excitation 
of all other meta-atoms via the scattered radiation. In 
this section, we combine the response of the meta-atoms 
to EM field and the expressions of the EM fields scat- 
tered by the meta-atoms in order to investigate how the 
radiation mediates interactions between the meta-atoms. 

We begin by examining the dynamics of a single driven 
meta-atom in Sec. IV Al There we show that when radia- 
tive losses are much weaker than the resonance frequency, 
a single meta-atom's dynamics reduce to those of the fa- 
miliar damped LC circuit in which the energy is lost to 
the scattered EM field. We then examine interactions 
between different meta-atoms in a collection of closely- 
spaced resonators. Due to the strong coupling between 
the EM fields and the current oscillations, the emitted ra- 
diation leads to the collective dynamics of the ensemble. 
In Sec. lVBl we explore the collective response of the sys- 
tem in the rotating wave approximation, in which each 
meta-atom's radiative emission rate is much less than its 
resonance frequency. We present an analysis for a more 
strongly interacting system outside the rotating wave ap- 
proximation in Appendix [C] 

In these treatments, we assume the spatial extent of 
each meta-atom is much smaller than the wavelength of 
EM field with which it interacts. As such, the radia- 
tion scattered from each meta-atom can often be approx- 
imated as that of electric and magnetic dipolcs oscillating 
in sync with one and other. For simplicity, when eval- 
uating the interactions between meta-atoms, we assume 
that the electric quadrupole and higher order multipolc 
contributions to the radiation of a single meta-atom are 
much weaker than the dipolc radiation and that they can 
be neglected. This is by no means a necessary approxi- 
mation. We could extend the general formalism to incor- 
porate multipolc-field radiation components in a multi- 
pole expansion. The dipole approximation, however, will 
provide an advantage in maintaining the tractability of 
the derivation of the collective mctamaterial response to 
EM fields. Moreover, in several practical situations, a 
unit-cell resonator of a metamaterial array may consist 
of two or more meta-atoms. Hence, in the dipole approx- 
imation to a single meta-atom, the unit-cell resonator 
would still exhibit multipolc radiation contributions. The 
multipole fields radiated by unit-cell resonators are also 
weak in many cases. For instance, metamaterial samples 
consisting of asymmetric split ring metamolccules have 
been experimentally employed in the studies of collective 
resonator response i 32 ' 34 i 35 ' 52 In an asymmetric split ring 



metamolecule the generated quadrupole field is notably 
suppressed when compared to the corresponding dipolar 
fields 

The electric and magnetic dipole moments produced 
by the current oscillation in meta-atom j are 

dj = Qjhjdj and m 7 = IjAjAj , (77) 

respectively. These are given in terms of the charge Qj 
and the current Ij of the meta-atom. The geometry- 
dependent proportionality coefficients hj and Aj have 
units of length and area and are defined such that 

hjdj = j d 3 r Pj(r) and Ajihj = J d 3 r Wj(r). (78) 

The unit vectors dj indicate the orientation of the electric 
dipole while the unit vectors rhj indicate the orientations 
of the magnetic dipoles. The distributions Pj(r) and 
Wj(r) [see Eq. ((9|)] represent the spatial profile of the po- 
larization and magnetization densities in terms of Qj and 
Ij. While, generally, the current resulting from the polar- 
ization density [the first term in Eq. (|10bj) ] contributes to 
the magnetic dipole, the polarization and magnetization 
densities (and hence the mode functions) that lead to a 
given charge and current distribution are not unique 42 
We have therefore chosen for each meta-atom j, p_j, wj 
and the position vector Rj such that the contribution of 
the polarization current to the magnetic dipole moment 
about Rj is zero. 

To facilitate an understanding of how the EM field 
influences the meta-atom dynamics, we consider a meta- 
atom's self-generated fields separately from the fields gen- 
erated externally. Consider the dynamics of a single 
meta-atom j interacting with the EM field. The meta- 
atom's equations of motion are given by Eq. Q38|). To 
isolate the dynamics arising from the self-generated field, 
we decompose the electric and magnetic fields into those 
generated by meta-atom j - Eg,j and Bgj - and those 
generated externally to meta-atom j, Ej. oxt and Bj. cxt . 
We then obtain the following relationship between the 
different contributions 

E A «t = E ta + Es J' ( 79a ) 

Bj.ext = B in + $^B S)J v. (79b) 

These external fields include contributions from the inci- 
dent field and the fields scattered by all the other meta- 
atoms in the system. 

In the previous section we derived the expressions for 
the scattered fields in terms of the polarization and the 
magnetization densities of the source medium. It was 
advantageous to represent the scattered fields in the fre- 
quency domain. We similarly analyze here the Fourier 
components of the dynamic variable Qj oscillating at fre- 
quency CI. As we did with the emitted fields, we find it 
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convenient to decompose the meta-atom variables 



Qi(*) = o5 +) (*)+or(*)> 



(-)/ 



^(t) = 4 + \t) + 



4~ ] (t), 



(80) 
(81) 



into their positive and negative frequency components, 
with 



^ _) (*) = 



(82) 



The positive and negative frequency components for 
these variables are defined such that, for a given fre- 
quency f2 



Qf\ft) = Qj(ft)Q(±ft) 



</>)=> (n) = </>j(Q)e{±n). 



(83a) 
(83b) 



With frequency components of Qj and <f>j defined in this 
way, the positive (negative) frequency components of the 
dynamic variables are driven exclusively by the positive 
(negative) frequency components of the EM fields. Since 
the metamatcrial system we consider in this model is 
linear, the equations of motion in Fourier space become 
the algebraic relationships between Fourier components 
of a common frequency fi, 



<j)j(ft) - $j, se lf(») - $j,cxt(») 
h 



-ift(f>j(ft) = £j j8e lf(fi) +£j,ext(fy, 



(84a) 
(84b) 



where £j, se if and se if are the self-generated EMF and 
flux, respectively, while £j,ext and 3?j iex t are the EMF and 
flux generated externally to meta-atom j. The current re- 
lates to conjugate momentum and magnetic flux through 
Eq. (fT9| , and the equation of motion for Qj , Eq. (|84ap , is 
nothing more that the statement that the rate of change 
of Qj is the current Ij . This translates to the relationship 
between frequency components — iftQj(ft) = Jj (f2) . The 
EMF and magnetic flux contain the external driving in- 
duced by the external EM fields as well as driving induced 
by the field that the current oscillation itself generates. 

The externally applied EMF and magnetic flux are 
given explicitly in terms of the externally generated fields 
as 



£ jt e>a.(ft) = / d 3 rPj{r) ■ Ej,ext(r,n), 



$J,ext(fi) 



(85) 



J d 3 rw,(r)-B JiCXt (r,fi). (86) 

When the external fields vary slowly over the volume of 
meta-atom j, £j. cxt and ^ext reduce to a direct driving 
of the meta-atoms' electric and magnetic dipolcs, respec- 
tively 



£j,ext (ft) 



hjdj 



E 



j,ext 



(Rj,n), 



j,ext 



Ajrhj ■ Bj iCXt (Rj, ft), 



(87a) 
(87b) 

where Rj is the position of the meta-atom. The external 
EMF and flux mediate the interactions between distinct 
meta-atoms which we will discuss in Subsection IV Bl and 
Appendix [CJ 



A single meta-atom interacting with the EM 
field 



Before investigating how scattered EM fields facilitate 
interactions between meta-atoms, we first shed light on 
how the meta-atom's field influences the evolution of the 
meta-atom itself. This is done by studying a single, 
isolated externally driven meta-atom. We will present 
expressions for the self-generated fields' contribution to 
both the EMF and the flux. When the spatial extent 
of the meta-atom is much less than a wavelength, the 
self-induced EMF can be written in terms of an effective 
self-capacitance, and the magnetic flux can be written in 
terms of a magnetic self-inductance. We thus show how 
each meta-atom can be treated as a radiatively damped 
LC circuit which is driven by external fields. This anal- 
ogy allows us to define slowly varying normal variables 
and derive their dynamics. 



Self-induced EMF and magnetic flux 



The EMF and the magnetic flux represent reactions of 
a meta-atom to EM fields generated by the meta-atom 
itself, as well as to external fields. Self-generated electric 
and magnetic fields provide a major contribution to the 
EMF and magnetic flux, respectively. We define the self- 
generated EMF and flux as 



£j, so if(ft) = J d 3 r j Pj{vj)-'E Sij {r j ,n), 
$ i>se if(0) = j d 3 rw,(r) • B SlJ -(r,n) . (89) 



The self-generated fields of meta-atom j, i.e., the fields 
Esj and Hsj = Bs,j//xo — Mj scattered from meta- 
atom j, at a frequency f2 are given in [Eq. (|63|) ]. From 
the expression for Eg j [Eq. (|63a|) ] . we obtain the self- 
induced EMF [Eq. (|88[) ] in terms of the radiation kernels 
[Eqs. dUl) and ([67j)] 



£j, self 



d 3 r / d 3 r' 



47TC0 

p J (r)-G(r-r',0)-p 3 (r')g j (0) 



1 



+ -p j (r)-G x (r-r',ft)-v, j (r')I(ft) 



(90) 



Similarly, the self-generated flux is obtained from the ex- 
pression for Hgj [Eq. (|63b[) ]. and is given in terms of the 
radiation kernels as 
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2. Self- capacitance and self -inductance 



$3 )Se lf(n) 



47T 



d 3 r |w,(r)| 5 



dV 



w i (r)-G(r-r',n).w j (r')/ i (n) 



-cw J -(r).G x (r-r',n).p j (r / )Q i (f2) (91) 



The first term of Eq. ([91]) arises because the flux is de- 
fined in terms of B = /i (H + M) rather than H, whose 
scattered field components are determined by the radi- 
ation kernels. This results in different contact terms in 
Eqs. flU and ([9ljl. 

Because we have assumed that the meta-atoms are 
much smaller than the wavelength, we may expand the 
radiation kernels to lowest order in kr and thus approx- 
imate the self-interactions in the near field limit. Since, 
in this limit, G x (r, fi)/G(r, f2) ~ kr -C 1, we neglect the 
contribution of G x to the self-interaction. To leading 
order in kr, we have the positive and negative frequency 
components of the radiation kernels 



RcG (±) (r,ft) 



3f f 



ImG (±) (r,!l) w ± 



2 



47T 
3P 



(92) 
(93) 



The long wavelength approximation allows us to sim- 
plify the expressions for the self-generated EMF and flux 
[Eq. ([90)) and ([91]) ] by neglecting the contributions of 
G x - This approximation implies that the self- induced 
EMF is directly proportional to the charge Qj , and that 
the self-induced magnetic flux is directly proportional to 
the current Ij . We can thus draw an analogy between a 
typical meta-atom and a standard LC circuit where the 
charge Qj and current Ij are related to £j ySC \f and "Itself 
through an effective capacitance Cj and magnetic self- 
inductance Lj- From Eqs. ([92]) and ([93]) . the positive 
and negative frequency components of the EMF and flux 
arising from the meta-atom's self-generated field become 



i 



C 



L< M >±i 



67Te J 

6w 



(n), 



(94) 
(95) 



In addition to the capacitance and inductance the EMF 
and flux have respective imaginary contributions that, 
as we shall see later, represent dissipation of the current 
oscillation due to radiation being emitted away from the 
meta-atom. The self-capacitance Cj is given by 



1 [ d 3 r \PMl-^ / d V /^3(p 3 (r).n)(n.p,(Q) p j( rO. Pj (r) 
Cj J 3e 47re i J |r-r'| 3 

with n = (r — r') / |r — r'|, and the magnetic self-inductance is 

,(m)_2mo [j3 mlm . M Pj_»o f ^3^/ 3 ■ n) (n ■ w 3 (r')) - w,(r') ■ w,-(r) 

|r — r' 



L) m > = Z?L I d A r |w,(r)| z + ^ / d A r / rfV y Jy ' ,y , ,,3 — J , (97) 



r 



In essence, excitation of the dynamic variable Qj pro- 
duces a distribution of electric dipoles (polarization den- 
sity) proportional to the mode function Pj(r). In 
the long-wavelength approximation, this distribution of 
dipoles produces a quasi-static electric field in the vicin- 
ity of the meta-atom generated by the real part of the 
radiation kernel ReG [Eq. ([92]) ] . The current oscilla- 
tion interacts with itself via the near field electric dipolc- 
dipole interactions, resulting in the effective capacitance 
C 3 appearing in the self- induced EMF [Eq. ([94]) ]. Simi- 
larly, a nonzero current Ij produces a distribut ion of mag- 
netic dipoles (magnetization density) proportional to the 
mode function Wj . The current oscillation then interacts 
with itself via the near field magnetic dipole-dipole inter- 
actions, resulting in the magnetic self-inductance ij M ^ 



appearing in the self-induced flux [Eq. ([95]) ]. 

Because the self-induced flux [Eq. ([95]) ] is proportional 
to Ij , we find it convenient to express the conjugate mo- 
mentum = Ijlj + <f>j [Eq. (|19[) ] in terms of a total 
sclf-inductance 

Lj^lj + L^. (98) 

This self-inductance includes contributions from both the 
magnetic and the kinetic inductances. When we include 
contributions from both the self-generated flux [Eq. (|95jl ] 
and the external flux [Eq. (|86p], the conjugate momen- 
tum for meta-atom j is given in terms of the total self- 
inductance by 

# ] (0) = (lj ± *^f^J if + *gL(n). (99) 
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This relation will be useful in determining the meta-atom 
equations of motion. 



3. Equations of motion for a meta-atom interacting with 
its self- generated fields 

Having determined how the self-scattered fields affect 
the EMF and flux, we now determine a closed set of equa- 
tions of motion for the meta-atom's dynamic variable and 
conjugate momentum. The rate of change of the dynamic 
variable Qj is given by the current Ij . Solving Eq. (|99|) for 
thus allows us to determine an equation of motion 
for Qj in terms of its conjugate momentum and magnetic 
flux generated by the external field. Further, substituting 
the EMF from Eq. ([94]) into Eq. ([84]) provides the cor- 
responding equation of motion for <pj. Explicitly these 
equations of motion are given in the frequency domain 
as 



(±) 



(±) 



1 



ck\ 



M,j 



Uj Df> 



If 



4> 



(±) 

i,ext 



(100) 



**j ' °j,ext 



where, as we demonstrate later, 

1 



(101) 
(102) 



is the single meta-atom resonance frequency, k = |0|/c 
is the wavenumber of the field frequency component, 



67T60C 3 

is the emission rate due to electric dipole radiation, 



(103) 



67TC 3 ii 



(104) 



is the emission rate due to magnetic dipole radiation, and 
Df\il) = l±i(±Y^* (105) 



arises from the inversion of Eq. ([99]) . The interaction 
of the meta-atom with its external fields are parameter- 
ized by hj and Aj [Eq. ([78]) ] and hence by the radia- 
tive emission rates Tej and Tmj- This is made clear in 
the point dipole approximation where we have the exter- 
nal EMF and magnetic flux which drive the meta-atom 
[Eq. ([87])]. From Eqs. (fT03]) and ([T04]) . one can infer that, 
when the meta-atom geometry is altered such that the 
self-capacitance and self- inductance remain constant, an 
increased interaction strength of the meta-atom with the 
external field corresponds to increased radiative emission 
rates. 



4- A meta-atom as an LC circuit 

If we neglect the radiative damping and consider a 
meta-atom interacting exclusively with its self-generated 
field, its dynamics are nothing more than those of an LC 
circuit with resonance frequency luj, which in the time 
domain satisfies the equations of motion 



d_ 

dt 



Qj(t) 



L 



The meta-atom normal mode variables 




(106) 



(107) 



and (3* evolve with eigenfrequencies ujj and —ujj, respec 
tively 



/3j(t) = exp(-iw J -t) / 9 i (0). 



(108) 



The collective dynamics within the metamaterial, of 
course, arise from the interaction of each meta-atom with 
its external field, necessitating the inclusion of radiative 
losses Tej, Tmj- But, as we will see later in this section, 
the presence of radiative interactions not only results in 
energy being carried away from the meta-atom by the 
radiated field, but also allows the meta-atom to be driven 
by fields scattered from other meta-atoms. 



5. The meta-atom normal oscillator variables 

The variables (3j represent cigenmodes of a single meta- 
atom in the absence of interactions with the external 
fields. The presence of these interactions perturbs the 
single meta-atom dynamics. Since the incident EM field 
driving the metamaterial oscillates at a central frequency 
f2o, it is convenient to analyze the effects of these pertur- 
bations using the slowly varying normal oscillator vari- 
ables 



b 3 (t) 




(109) 



The oscillator variables satisfy the Poisson brackets 





{b(t),b(t)} 

{&,•(*),&;,(*)} 



{b*(t),b*,(t)} 
-i5j,j>. 



(110a) 
(110b) 



One can recover Qj and 4>j by solving the system of equa- 
tions formed by Eq. (|109|) and its complex conjugate. 
This yields 



>(*) 



1 

V2 



(e- inot b 3 (t) + e mot b*{t)) 



= -i-^{e- iQ ° t b j (t)-e in ° t b*(t)) 



(111) 
(112) 
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As the incident electric field may consist of a range of 
frequencies around 51o reflecting its variation in time, it 
is necessary, in general, to examine the frequency com- 
ponents of the oscillator variables and how they are re- 
lated to those of Qj and <pj. The Fourier components, 
for > 0, of Qj and 4>j are given in terms of the normal 
variables as 



(Q) 



1 



(b J (5) + b*(-S-2n )), 



y/ U 3 L 3 



-i—( bj (6)-b*(-6-2n )) 



where 



(113) 
(114) 

(115) 



The negative frequency components of Qj and <j)j [given 
in terms of their positive frequency components in the 
time domain in Eq. (1821)]. when O < 0, can be ob- 



tained from the relations, Q^ '(0) 



and 



*«-'(o)= *; +, (-si) 



3« + », 



The RWA essentially assumes that all the dynamics 
are dominated by the frequency £1 . Because the RWA 
implies <5f2, \uij — f2o| <C Qo, we can approximate the 
quantities (Q/ojj) 3 appearing in the equations of motion 
[Eqs. ([TOO)) and (TlOTj) ] as (fl/(Jj) 3 « 1. In these limits, 
the equations of motion for the frequency components of 
Qj and fa [Eqs. (TT0T)|) and ([TOT]) ] yield the relationship 
for the normal variables 



i<%(<5) = 



—i (ujj — Qq) — 



(118) 



where the detuning of the meta-atom resonance ujj from 
the frequency of the driving field Oo manifests itself as an 
oscillation of the normal variable bj at frequency uij — Slo, 
while electric and magnetic dipole radiation emanating 
from the meta-atom results in the damping of bj at a 



rate Te. 



The forcing function combines driving 



of the current oscillation by the external electric field via 
the EMF and the external magnetic field via the flux, 
and is given by 



iuj$ jt ex t (5) 



(119) 



6. Dynamics in the rotating wave approximation 

Radiative damping and driving of the meta-atom by 
external fields alter the current oscillation represented 
by the normal variable bj. The interactions leading to 
these effects are often sufficiently weak that we can re- 
gard their influence as a small perturbation. We con- 
sider this weak interaction limit here and in Sec. IV Bl 
where we examine the collective behavior of the meta- 
atoms comprising a metamaterial. We thus assume that 
bj varies slowly with respect to the dominant frequency 
fio and neglect the fast oscillating components, i.e., we 
set bj(—6 — 20o) = for \S\ <C ^o- The mode vari- 
ables bj are then proportional to the slowly varying en- 
velope of the positive frequency components of the dy- 



namic variables Q^ and their conjugate momenta r 
Neglecting fast oscillating components of bj is known as 
the rotating wave approximation (RWA), and is valid in 
the limit T^.j, Tmj", |^o — ujj\,SQ <C f^o, where 60, <C Oo 
indicates a narrow bandwidth of the incident field. 

In the RWA, the meta-atom driving forces, i.e., the 
EMF and flux, can be expressed in terms of their slowly 
varying envelopes £j and $j (t) defined such that 



(+) 



^/LJjLj 

$( +) (t) 



e- iQot Sj{t) 



(116) 
(117) 



where the overall factor of ^uijLj was included for con- 
venience. 



7. The meta-atom as a driven, RLC circuit 

Here, we show that in the RWA, a meta-atom behaves 
as a damped, driven RLC circuit interacting with the 
external driving field. A source of loss that is typically 
present in a meta-atom which we have thus far neglected 
is the ohmic losses due to resistance to current flow within 
the meta-atom. We include the effects of this resistance 
phenomenologically through the addition of the ohmic 
loss rate ^o,j to the radiative damping rate. The ap- 
proximations leading to Eq. (j!21[) are still valid provided 
that Tq Hq. The total meta-atom damping becomes 



Tj=T 



r 



Mj 



r 



o,j- 



(120) 



We obtain equation of motion for bj in the time domain 



from Eq. (|118p by multiplying by 



-i5i 



and integrating 



over the bandwidth of the external field — SQ < 5 < 5Q. 



dbj_ 
dt 



-i (ujj - O ) 



bj{t)+f : 



j'.ext 



(t) (121) 



When the incident field is of finite duration, i.e., 
E in (r, ±oo) = B in (r, ±oo) = 0, bj satisfies Eq. (|12ip with 
the initial condition bj{—oo) = 0. 

The interaction of the meta-atom's current oscillation 
with its self-generated EM fields cause the current mode 
to oscillate at the resonance frequency ojj [Eq. (|102[) ] 
analogous to that of an LC circuit. When a meta-atom 
current oscillation produces net electric and magnetic 
dipole moments, this oscillation can be driven by exter- 
nal fields as manifested by the term /j, e xt(i) m Eq. (|12ip . 
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Radiative and ohmic losses act as a resistance within the 
meta-atom, and the external EMF £j, C xt and <frj. cyc t pro- 
vide the driving. 

The dynamics of this effective RLC circuit can be de- 
rived from the effective Hamiltonian 



HeS,j = H^Cj + Hdamp,j + Kxt ,j 



(122) 



where -Hlcj is the effective Hamiltonian for an un- 
damped LC circuit, 



(123) 



the damping is provided by adding an imaginary term to 
the effective Hamiltonian 



-iTbjbj, 



(124) 



and the interaction with the external field is provided by 

C.c. , (125) 



Vextj — °j,cxt^j j,ext 



Li 



The physical significance of the interaction term becomes 
clearer in the dipole approximation. When we neglect the 
spatial extent of the meta-atoms, the interaction poten- 
tial with the external fields beomes 



cxt,j 



(+) 

j>ext 

(+) 
j,ext 



(Rj,i)-mJ _) (i) + C.c. , (126) 



where d^\t) = hjQj (t) is the electric dipole of the 
meta-atom, and 



A 
L 



(127) 



is an effective magnetic dipole of the meta-atom. To 
understand why m'j can be interpreted in this way, 
consider the conjugate momentum expressed in terms of 
the self-inductance [Eq. ([99]) ] in the limits of the RWA 
(namely CI/uj ss 1) 



Li 1 ± 



r 



(±) 



(±) 

j',ext 



(128) 



Because Fmj "C Wj, when the self-induced magnetic flux 
dominates that generated by external fields, the conju- 
gate momentum is related to the current by 



(±) 



(129) 



and vaj sa nTj^ is approximately the magnetic dipole 
created by the current oscillation in meta-atom j. The 
effective interaction Hamiltonian [Eq. (j!26[) ] accounts for 
the energy of the meta-atom electric dipole interacting 
with externally generated electric fields and the meta- 
atom's magnetic dipole interacting with externally gen- 
erated magnetic fields. 



The energy lost due to radiative damping is carried off 
by the scattered fields. The external fields contributing to 
the interaction V e xt,j include fields scattered from other 
meta-atoms in the system. In the following subsection we 
will explore how these scattered fields drive and influence 
the dynamics of the meta-atoms. 



B. Collective interactions in the rotating wave 
approximation 

In this subsection, we examine in detail how the fields 
emitted externally to meta-atom j drive the excitation in 
that meta-atom. In particular, we will see how the fields 
emitted or scattered from the ensemble of meta-atoms 
mediate interactions between them. The EM field gener- 
ated externally to each meta-atom has two components: 
the incident field, and the fields scattered from all other 
meta-atoms in the system. The incident field impinges 
on the metamatcrial driving all of its constituent meta- 
atoms. Each excited meta-atom, in turn, radiates an EM 
field which can drive other meta-atoms while undergoing 
multiple scattering between different resonators. In or- 
der to calculate the response of the metamaterial array 
to incident EM fields, we need to consider these multiple 
scattering processes, which produce a coupling between 
meta-atom current oscillations. For near-resonant fields, 
recurrent scattering events in which the field scatters off 
the same meta-atom multiple times dramatically affect 
the potentially strong coupling between closely-spaced 
resonators. 

Here we will derive a coupled set of equations for the 
meta-atoms where all the multiple scattering processes 
are fully incorporated in the EM field induced interac- 
tions between the meta-atoms. We will then examine 
how the coupling can lead to a cooperative response of 
the metamaterial to the incident field via excitation of 
collective modes of current oscillation. Such modes can 
have either superradiant character, where the interac- 
tions enhance the radiation emitted from metamaterial, 
or a subradiant character, where the radiation remains 
trapped as it repeatedly scatters between meta-atoms 
leading to a suppressed collective radiative emission rate. 

In order to derive a coupled set of equations for the 
meta-atoms where the interactions are mediated by the 
EM fields we consider the meta-atom mode variables 
bj and investigate their dynamics within the RWA. As 
stated in Sec. IV Al in order for the RWA to be valid, we 
assume that the emission rates satisfy T-gj, Tyi j, To.j <C 
fio and that the driving field's bandwidth and its detun- 
ing from meta-atom resonance are small compared to the 
frequency of the driving field, i.e., <5f2, — flo\ <C f^o f° r 
all meta-atoms j. In these limits, the external field inter- 
actions act as a slow perturbation on the fast oscillations 
caused by the meta-atoms' self-generated fields. 

A meta-atom j experiences driving from the external 
electric and magnetic fields. These fields induce EMFs 
and fluxes, which by Eq. ()121[) . impact the dynamics of 
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the current oscillation. The driving originates from both 
the incident field, and from the fields scattered from all 
other meta-atoms f 7^ j in the system. As such, we 
decompose the EMF and flux into those generated di- 
rectly by the incident driving, £j; m and "i^in, and those 
induced by fields arriving from meta-atom j', £j j> and 
&ji>- Explicitly, 



(130) 
(131) 



The incident field directly drives each meta-atom, induc- 
ing a forcing term 



1 

3 



(132) 



while the scattered fields produce a coupling between the 
resonators. Below, we will show that in the RWA, the 
scattered fields emanating from meta-atom j 1 are propor- 
tional to the amplitude by of the oscillation in meta-atom 
j', and therefore that Ej^y and $j- 3 -< are proportional to 
by. We will find that, by virtue of the scattered fields, 
the dynamics of the individual meta-atoms are coupled. 
The ensemble will exhibit collective modes of oscillation, 
each with its own frequency and radiative decay rate. 

Because the incident field has a narrow bandwidth 
around a frequency Oo, we find it convenient to define 
slowly varying quantities to describe the dynamics of 
the system. For any vector field F(r, i) = F^ + ^(r, t) + 
F^'(r,t) with positive and negative frequency compo- 
nents F^" 1 ") and F^ - ), respectively, unless otherwise spec- 
ified, we define the slowly varying envelope F(r, t) of the 
field such that the positive frequency component 



pW(r,t) = e- tnot F(r,t) , 
or equivalently in frequency space 

F<+>(r,n)=F(r,*) , 



(133) 



(134) 



where again S = O — O [Eq. (|115|) ]]. For the charge 
and conjugate momentum on meta-atom j, we define the 
scaled slowly varying quantities Qj and <pj such that 



Q\ + \t) 



y*UJyCj 



e~ inot Qy{t) 



- inot 4>yit) 



(135) 
(136) 



In the RWA Qj and <j>j arc trivially related to the normal 
variables by 



/(*) 



V2 



(137a) 
(137b) 



Outside the RWA, bj contains fast oscillating compo- 
nents whose origins we discuss in Appendix [C] In this 
subsection, however, we will assume that Eq. (|137p holds. 
We also define the scaled current such that 



/^/(+)( t )= e -^ J-(t) 

,0j 



(138) 



The relative scale factor of the current was chosen so 
that, for a frequency 5 <§C Slo, the Fourier components of 
cf)j and Ij are related by 

4>y(5) = Dj(Qo + 5)Ij(S) + $ jjex t(5) • (139) 

The quantity Dj [Eq. (|105[) ] serves as the dimension- 
less complex self-inductance. Because we have assumed 
Tmj- "C fio in the RWA, the quantity n 3 - « 1. 

Next we will determine the contribution of the fields 
scattered from each meta-atom j 1 to the normalized 
EMF, Ejy, and flux, &j t y, of meta-atom j. We express 
the scattered fields from the meta-atom j' in terms of the 
normalized variables Qy and Ij' . We assume the band- 
width of the incident field is sufficiently small that the 
time scale over which the fields vary, 1 / SO,, is much longer 
than the time it takes for light to propagate across the 
mctamatcrial sample. We then obtain the slowly varying 
scattered fields by substituting D,q for n in the radiation 
kernels, G and G x [Eqs. (|63a|) and (|63b|) ]. and exploit 
Eq. (fT54]) to obtain 



E 



s,j' 




VTijQj, J d 3 ry G(i 



.fin 



+ V r M,y Iy I d 3 ry G x (r - vy , fi ) 
and 



Pi' fa') 



w i'( r i') 

Ay 



(140) 



H 



S,j' 



3 j4_ / no 

2y 6iT[io \ w i' 



3/2 



^J'iy J d\y G W (r - Vy , flo) • 



v/TVQ,-, / d 3 r 3 , G ( x +) (r - iy , n ) 



Pj'( r j') 



(141) 

The amplitude of the electric and magnetic fields emitted 
by the electric dipole of meta-atom j' , driven by Qj', 
scale with •v/Tej/. Similarly, the fields emitted by the 
magnetic dipole of meta-atom j' , driven by ly , scale with 

These scattered fields provide a portion of the slowly 
varying EMF, 

1 



,/ ;! r ; p ; ir/i- Ks.^r,; (142) 
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and flux 



Mo 



d 3 r,w,(r,).H Sj ,(r,) (143) 



at meta-atom j. Substituting Eqs. (jMOj) and (|141|) into 
the expressions for EMF and flux gives 



°3,3 



VT^r Ed , [Q E (n )] jtj , o y 



where the matrices Ge, Sm, and G x determine how the 
meta-atoms' geometries and relative orientations influ- 
ence the respective contributions of the scattered electric 
fields to the EMFs, the scattered magnetic fields to the 
fluxes, and the scattered electric (magnetic) fields the 
fluxes (EMFs). These matrices have zero diagonal ele- 
ments and off diagonal elements given by 



$.,=1 

^3,3 



Q 



^3 KV^W' 



x/FmjFmj-' [Qm{P,q)]-_-, Ij> 

(145) 



V^VT^W [Gxfoo^jj' Qj' 



[GE(Q)] jtJ , 



G(r 3 --r 3 v,n) 



Pj'( r j') 



hy 



™3'( T 3') 



Ay 



G x (r 3 - ry,Ct) 



W 3'( r 3') 



(146a) 
(146b) 
(146c) 



r 



When the separation between two meta-atoms is much 
greater than the spatial extent of the individual elements, 
these geometrical factors depend exclusively on the rela- 
tive positions and orientations of the meta-atoms' electric 
and magnetic dipoles. Explicitly, in that limit, 



[9E(n)] jd > 

[Ou(myy 
[<3xm jt y 



¥<■ 

3 „ 
-m, 

2 3 



G(R 3 --Rj/,n).dy (147a) 
• G(Rj - Ry , fi) • uiy (147b) 
G x (Rj - Ry , 0) • ihy ■ (147c) 



The contribution of the electric field scattered by meta- 
atom j' to the EMF, £j y , scales with the geometric mean 
of the electric dipole emission rates of the two meta- 
Similarly, the magnetic field of cl- 

with a strength 
When the meta-atoms are 



atoms, wTejTeJ' 



cment j' contributes to the flux §>jj 
proportional to WT^jT^j 
sufficiently far away from one and other, the electric field 
emitted by the magnetic dipoles and the magnetic field 
emitted by the electric dipoles provide a significant con- 
tribution to Eyy and ^yy that scale with a/FejTmj' 
and ^/FmjTejv, respectively. 

We have set out to obtain coupled equations of motion 
for the meta-atom normal variables bj mediated by the 



EM field. We have obtained contributions to the EMF 
and flux that are driven by charges Qj and currents Ij. 
However, only Qj and conjugate momenta (f>j are triv- 
ially related to these normal variables [Eq. (|137j) ]. The 
current, on the other hand obeys the more complex rela- 
tionship 



h> = 



V2 



(148) 



One can thus use Eqs. (|137|) and (|148p to express £j t y 
and in terms of the normal variables by. We note, 



j'- 

contains contribu- 

~^/ujj 



however, from Eq. (|145|) . that §j> 

tions that scale as ^TMj'^M.j" I and ^TM.j'^E.j" / j 
which under the conditions of the RWA, are much less 
than 1. Furthermore, Ij contains a contribution from the 
incident field flux $j,i n - The contribution of the incident 
flux to Ij can also be ignored to lowest order since it is 
about Tj/ujj times the direct contribution of the incident 
flux to the direct driving fj; m - So, to determine Sjji and 
to lowest order in Tj/uij, we therefore exploit the 
approximate relationship Iy rj §y w —iby /V2. 

Having computed the contributions of the scattered 
fields to the EMF and flux of an individual meta-atom, 
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wc find that these scattered fields produce a coupling be- 
tween meta-atoms in the oscillator equations of motion. 
Substituting the EMF and flux into Eq. (fT2"Tj) . wc find 
the evolution of the column vector b of normal variables 
is governed by 



b = Cb + f ir 



(149) 



where we have introduced the following notation for b 
and for the driving fi n caused by the incident field: 



(150) 



The coupling matrix C is given to lowest order in 
Tej/loj', and T u ,j/oJj' by 
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Here the detunings of the incident field from the meta- 
atom resonances are contained in the diagonal matrix A 
with elements 



Moreover, the meta-atom emission rates are incorporated 
in the diagonal matrices Ye, Tm and Tq with elements 





= r Ej , 


(153a) 






(153b) 




= r o,j , 


(153c) 



respectively, and we have defined T = Ye + Tm + Yq- 

The interaction matrix C accounts for electric dipole- 
dipole interactions, magnetic dipole-dipole interactions, 
as well as interactions between electric and magnetic 
dipoles that arise from magnetic (electric) fields emitted 
by electric (magnetic) dipoles. The diagonal elements of 
C result from interactions with the self-generated fields 
and give rise to the meta-atoms' resonance frequencies 
and radiative emission rates. 

In the RWA, the dynamic equation [Eq. (|149|) ] encap- 
sulates all the multiple scattering processes between the 
different meta-atoms. These are described by the interac- 
tion terms in the matrix C, mediated by the scattered EM 
fields. The coupled set of equations implies a system of N 
meta-atoms possesses TV collective modes of excitation. 
These modes correspond to the eigenvectors of the ma- 
trix C. For each collective eigenmode we have collective 
radiative resonance linewidths and resonance frequencies 
that are represented by the eigenvalues of C. A strong 
coupling between the resonators can lead to a cooperative 



response of the mctamaterial sample to the EM fields, re- 
sulting in collective decay rates which are substantially 
different from those of a single, isolated meta-atom. The 
interactions can either enhance radiative emission, pro- 
ducing a superradiant mode, or suppress emission, yield- 
ing a subradiant decay rate. We will illustrate the effect 
of a cooperative response of a 2D metamaterial array in 
Sec. I VII by considering an example of closely-spaced split 
ring resonators. We find that even in a relatively small 
sample the strong coupling leads to a dramatic resonance 
linewidth narrowing of five orders of magnitude and to a 
broad distribution of radiative decay rates. 

In order to illustrate the coupling of an incoming field 
to collective modes, suppose the incident field is engi- 
neered so that it only excites the i th collective mode, and 
then is suddenly turned off. The collective excitation is 
then distributed over the sample according to the eigen- 
vector Vi of C. Due to the repeatedly scattered fields 
that couple the meta-atoms, the excitation oscillates at 
its resonance frequency given by the eigenvalue Aj, 

fii = fi - Im(Ai) , (154) 

and the amplitude of oscillations decay at a rate 

7< = -2Re(Ai). (155) 

as radiation leaks out of the collective excitation and en- 
ergy dissipates through ohmic losses. The vector of nor- 
mal variables then evolves as 



(152) b(i)ocexp{ -i(fii-O )- — tj 



(156) 



The nature of collective modes could also allow one 
to engineer a cooperative response of the metamaterial 
to the incident field, addressing linear combinations of 
modes by shaping the incident field's profile, or adjust- 
ing its frequency. Engineering of the collective response 
may then be used, for example, to excite isolated sub- 
wavelength hot spots in a metamaterial.— 



C. Concluding remarks 

In this section, we saw how the interaction of individ- 
ual meta-atoms with the EM field governs the collective 
dynamics of an ensemble of meta-atoms that make up a 
metamaterial. Each meta-atom experiences the influence 
of its current oscillation's self-generated field, the field 
incident on the metamaterial, and the fields scattered 
from all other meta-atoms in the system. We explored 
the influence of the self-generated fields in Sec. IV Al In 
the RWA, the self-generated field dominates meta-atom 
dynamics. Each meta-atom can be seen as an effective 
RLC circuit which experiences damping due to electric 
and magnetic dipole radiation carrying energy away from 
the meta-atom. On the other hand, fields generated ex- 
ternally to the meta-atom, i.e., the incident field and the 
fields radiated from all other meta-atoms in the metama- 
terial, drive the current oscillations in each meta-atom. 
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In Sec. lVBi wc saw how the fields scattered by each meta- 
atom mediate interactions between them. Fields emitted 
by one meta-atom drive the current oscillations in all the 
others, producing the dynamic inter-meta-atom coupling 
in Eq. t) 149[) . While Appendix [Cl develops a formalism 
to account for arbitrarily strong interactions, in this sec- 
tion we have gained a significant physical insight in the 
RWA in which we assume the mcta-atoms' interact much 
more strongly with their self-generated fields than with 
the fields generated externally. 

In the following section, we will apply this formalism 
to examine collective modes in an example metamaterial: 
an array of symmetric split ring resonators. This sys- 
tem will illustrate the vital role cooperative interactions 
can play in the dynamics of a metamaterial composed 
of closely spaced plasmonic resonators. A metamaterial 
of N resonators will have N collective modes of current 
oscillation, each with its own resonance frequency and ra- 
diative emission rate. Both of these quantities strongly 
influence how a given mode can be excited. The coopera- 
tive interactions lead to a broad distribution of collective 
decay rates indicating strongly superradiant or subradi- 
ant modes. 



VI. AN ENSEMBLE OF SYMMETRIC SPLIT 
RING RESONATORS 

In this section, we apply the formalism developed in 
this article to a metamaterial composed of split ring res- 
onators (SRRs). As the name suggests, these resonators 
are composed of loops with segments that have been re- 
moved. Owing to the curvature of the elements, current 
oscillations within SRRs can exhibit both an electric and 
a magnetic response. Variations of these resonators have 
been used to produce metamaterials which exhibit, e.g., 
negative indices of refraction*^ Here, we consider a par- 
ticular realization of the SRR in which a single ring is cut 
into two disconnected concentric circular arcs of equal 
length. We then study the SRR metamolecule by assum- 
ing that the halves each form a meta-atom that supports 
a single mode of current oscillation. The two halves could 
either oscillate in phase, producing a net electric dipole, 
or out of phase, producing a net magnetic dipole. 

In addition to active studies of metamaterial arrays of 
SRRs, there has also been an increasing interest in fabri- 
cating metamaterials consisting of split ring resonators in 
which the symmetry between the two disconnected halves 
has been broken, e.g., by making one of them longer. 
Sheets of asymmetric split ring resonators (ASRs) have 
been shown to exhibit transmission resonances^ corre- 
sponding to excitations in which all magnetic dipoles in 
the sheet oscillated in phase. The quality factor of this 
resonance, however, was shown to depend strongly on 
the number of ASRs in the system.— Furthermore, arti- 
ficially adjusted disorder in the positions of the unit-cell 
resonators was observed to destroy the resonance^ If in- 
teractions mediated by the EM fields were not important, 
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FIG. 1: (color online) A schematic illustration of a split ring 
resonator. An excitation in the left meta-atom (I) produces 
an oscillating electric dipole (indicated by the blue arrow) in 
the direction d and a magnetic dipole (indicated by the red 
arrow) in the direction — rh, while an excitation in the right 
meta-atom (r) produces an electric dipole in the direction d 
and a magnetic dipole in the direction rh. The meta-atoms, 
in the point source approximation, are separated by a vector 
u. When the meta-atoms are excited in phase, the electric 
dipoles reinforce each other and the magnetic dipoles cancel 
out. 

and the ASRs behaved independently, system size or po- 
sitional disorder of the system would have little effect on 
the metamaterial response to the EM fields. These ex- 
perimental observations provide ample evidence for the 
vital role collective interactions play in this particular 
metamaterial. 

Here we employ the formalism describing collective in- 
teractions to an ensemble of SRRs in the RWA. We de- 
scribe a single SRR in subsection I VI A[ while we examine 
the properties of collective modes of SRRs in a lattice in 
subsection IVI Bl 



A. The symmetric split ring resonator 

We begin by describing the interaction of a single SRR 
unit-cell resonator with incident EM fields. This partic- 
ular realization of an SRR metamolecule consists of two 
meta-atoms formed by two concentric circular arcs la- 
beled by j € {l,r} (for "left" and "right"), as shown in 
Fig. [TJ This metamolecule possesses reflection symmetry 
about a central plane. 

To illustrate this qualitative physical behavior of an 
SRR, we approximate the meta-atoms as two point 
sources separated by u = R r — R; (see Fig. Q]). The cur- 
rent oscillations in mcta-atoms produce electric dipoles 
with orientation d r = df = d associated with charge os- 
cillating between the ends of the arcs. Owing to the 
curvature of the meta-atoms, these currents also pro- 
duce magnetic dipoles with opposite orientations rh r = 
— ih; = ih. The generated electric dipoles lie in the 
plane of the SRR and are perpendicular to the displace- 
ment between the meta-atoms (d 1 u). The generated 
magnetic dipoles, on the other hand, point out of the 
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plane in which the SRR resides (rh _L u, d). Each meta- 
atom in isolation supports a single mode of oscillation 
with resonance frequency u>q. Here we consider a reso- 
nant driving with the frequency of the incident field sat- 
isfying fio = For simplicity, we also assume each 
element possesses identical radiative and thermal decay 
rates r E /M/o,z = r E /M/o,r = Te/m/o- 



In the RWA, the normal variables b r and bi [Eq. f| 109[) 
with j € {r, I}] describe the states of the right and left 
halves, respectively, of a single SRR mctamolecule in iso- 
lation. We may now apply the previously developed the- 
ory for the EM field mediated interactions between meta- 
atoms to a single SRR unit-cell resonator consisting of 
these two meta-atoms. According to Eq. (|149p . the nor- 
mal variables b r and 6/ are coupled by the EM fields as 



= C 



sun 



fr.i 
fl.i 



(157) 



Here Csrr denotes the specific coupling matrix in this 
case between the two meta-atoms, as described in detail 
below. The incident field impinging on the SRR pro- 
duces the driving terms fj t - m for each meta-atom j = l,r 
[Eq. (|132[) ]. Considering the meta-atoms as point emit- 
ters, the incident field excites their electric and magnetic 
dipoles resulting in the simplified driving terms 

, d-E(R n t) A m-B(Rr,t) 
} ri in =ih j=== uqA j=== — , (158) 
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-, d-E(R^) , m-B(Ri,t) 
//, in =ih j== h w A j== — • (159) 
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The quantity h is an effective length along which charge 
flows to form the meta-atoms' electric dipoles and is re- 
lated to Te through Eq. f| 103[) . Similarly, A is an effective 
area that indicates the strength of the magnetic dipolc 
interaction and is related to the magnetic dipole emis- 
sion rate Tm through Eq. 104[) : L is the self- inductance 
of each meta-atom. Once excited, each half of the SRR 
scatters both electric and magnetic fields. These fields 
then impact the other meta-atom, driving its electric and 
magnetic dipoles. Repeated absorption and re-emission 
of scattered fields produces a dynamic interaction be- 
tween the two halves of the SRR. From Eq. (|151|) . the 
coupling matrix governing the interaction is given by 

r ( -r/2 4 drG-rs\ n „> 

Csrr - ldTG _ fs _ r/2 J (160) 
where a single meta-atom has a total decay rate 



r = r E + r M + r Q 



(161) 



appearing in the diagonal elements of Csrr, and we have 
defined 



r = vTeFm, 
dr = r E - r M • 



(162) 
(163) 



Coupling between the two halves of the SRR, represented 
by the off diagonal elements of Csrr, arises from interac- 
tions between the meta-atoms' electric dipoles, the meta- 
atoms' magnetic dipoles, as well as a cross interaction 
between the electric dipole of one meta-atom and the 
magnetic dipole of the other. The strength of the electric 
and magnetic dipole-dipole interactions is proportional to 
the radiative decay rates Te and Tm , respectively. These 
dipole-dipole interactions also depend on the spacing be- 
tween the meta-atoms and the relative orientations of 
the dipoles. This geometrical dependence shows up in 
the factor 

G = ^d • G(u, fi ) • d = • G(u, fi ) ■ m. (164) 

Notice that because identical meta-atom excitations (i.e., 
when bt — b r ) produce parallel electric dipoles, but an- 
tiparallcl magnetic dipoles, the electric and magnetic 
dipole interactions work against each other; the strength 
of interaction arising from the geometrical factor G is pro- 
portional to Te — Tm- An additional interaction arises 
from the electric dipoles interacting with fields scattered 
from the magnetic dipoles and vice versa. The geometric 
mean of the radiative decay rates, T [Eq. (|162p j. governs 
the strength of this interaction. Relative orientations of 
the electric dipole of the left (right) meta-atom and the 
magnetic dipole of the right (left) meta-atom appear in 
the geometrical factor 



S= -d- G x (u, fi ) • m r 



(165) 



To analyze the collective modes of the SRR, we con- 
sider the dynamics of symmetric c+ and antisymmetric 
c_ modes of oscillation defined by 
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C± = -j= (b r ± bi) . 



(166) 



These symmetric and antisymmetric variables represent 
the cigenmodes of the SRR. From the dynamic equation 
[Eq. jH2])] and the SRR coupling matrix [Eq. (flBO]) ]. one 
finds 

j t c± = {--y T i A S rr) c± + F ± , (167) 
where an incident field produces the driving terms 
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F± = -j= (frM ± j),in) 



(168) 



The interaction between the elements produces the decay 
rates 7± and shifts the resonance frequencies of the sym- 
metric and antisymmetric modes by equal and opposite 
amounts, A SRR 



7± = T E (1 ± 2Im(G)) 

+ r M (lT2Im(G))- 
A S rr = -2Re(G) (T B - T M ) 



2f Re(5) + r 
- 2f Im(5) . 



(169) 
(170) 
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An analogy can be drawn between these mctamolcc- 
ular current oscillations and atomic or molecular en- 
ergy levels . 48 ! 49 The symmetric and antisymmetric modes 
have respective resonance frequencies uiq + Asrr and 
ujq — Asrr. When excited, the symmetric mode decays 
at a rate 7+, while an excitation of the antisymmetric 
mode decays at rate 7_ . 

Excitation of the symmetric mode (c + ) produces a net 
electric dipole since the individual meta-atom electric 
dipoles oscillate in phase while the meta-atom magnetic 
dipoles approximately cancel each other out. Similarly, 
excitation of the antisymmetric mode (c-) produces a 
net magnetic dipole and the net effect of the electric 
dipole approximately cancels out. The symmetric and 
antisymmetric excitations will thus be referred to elec- 
tric and magnetic dipole excitations, respectively. When 
the spacing between the arcs u <C A, the decay rates 
simplify to 



7+ 
7- 



2r E - 
2r M 



-r 



(171a) 
(171b) 



The electric mode loses energy via electric dipole radia- 
tion, while the magnetic mode emits magnetic dipole ra- 
diation. In the absence of magnetic dipole interactions, 
the symmetric and antisymmetric modes are analogous 
to superradiant and subradiant states in a pair of closely 
spaced two-level atoms: when the two-level atoms are 
excited in phase, the radiative emission rate is enhanced, 
and it is suppressed when the atoms are excited out of 
phase. Furthermore, in the SRR metamolecule the elec- 
tric and magnetic modes are driven purely by the electric 
and magnetic fields, respectively, with F+ cx d ■ Ei n (R, t) 
and F_ oc rh • Bi n (R, t), where R denotes the center of 
mass coordinate of the SRR. 

When more than one SRR is present, radiation emitted 
from one SRR impacts and drives oscillations in another. 
The resulting interactions produce collective modes of os- 
cillation for the whole system. We examine this collective 
behavior in the following subsection. 



B. Collective modes in an ensemble of symmetric 
split rings 

Having discussed how EM field induced interactions 
arise between two meta-atoms in a single SRR meta- 
molecule, we now explore how a collection of meta- 
molecules can behave in concert when brought together 
to form a metamatcrial. As an example we consider a 
2D N x x N y array of SRRs arranged in a square lattice 
with lattice vectors ai = ae x and a.2 = ae y . This finite 
array resides in a region with free space (as opposed to 
e.g. periodic) boundary conditions. A single SRR oc- 
cupies each unit cell of the lattice. They are oriented 
such that symmetric oscillations produce electric dipoles 
along the direction d = e^, and antisymmetric oscilla- 
tions produce magnetic dipoles pointing out of the lattice 
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FIG. 2: (color online) The distributions of collective radiative 
emission rates in a 33 x 33 square array of SRRs. The distribu- 
tion is represented as a histogram of the log 10 of the collective 
emission rates. Panel (a) shows the distribution for an ensem- 
ble of SRRs with lattice spacing a = 0.5A. Panel (b) shows 
the distribution when the lattice spacing a — 1.4A. When 
the lattice spacing is larger, interactions become weaker, and 
cooperative effects are diminished, leading to a narrower dis- 
tribution of collective mode line widths. 



in the direction m = e z . In this section, we quantify the 
collective interactions by examining the collective eigen- 
modes of the system and showing how the interactions 
can lead to strongly modified radiative emission rates. 
We also illustrate from this model how a subwavclcngth 
inter-molecular spacing enhances the collective behavior 
of the system. In particular, we find that a subwavclcngth 
lattice spacing produces a much broader distribution of 
subradiant and superradiant collective decay rates. 

While an SRR in isolation possesses two modes with 
two collective resonance frequencies and two decay rates, 
the presence of interactions in an ensemble can produce a 
broad distribution of collective linewidths. The lattice of 
N x x N y SRRs possesses 2N x N y collective modes of oscil- 
lation, where the i th mode corresponds to an eigenvector 
Vj of the interaction matrix C [Eq. (|151|) ] . The resonance 
frequency of this collective mode is shifted from Hq by 
Si = Oj — £l and has a collective decay rate 7*. These 
are given in terms of the mode's eigenvalue A^ as 



Si = -ImAj , 
ji = -2ReA i; 



(172a) 
(172b) 



respectively. Here we consider an ensemble of SRRs 
whose elements have equal single-meta-atom electric and 
magnetic decay rates = Tm, and we take the sep- 
aration between constituent meta-atoms of an SRR to 
be u = 0.12A. Because the thermal losses are equal in 
all meta-atoms, their presence would add to the decay 
rates of each collective mode equally. Since here we are 
interested in how interactions modify collective radiative 
decay rates, we take the ohmic loss rate to be zero in this 
section. 

We numerically calculate all the eigenmodes of the sys- 
tem that are modified by the multiple scattering pro- 
cesses. Figure [5] illustrates how interactions mediated by 
the EM field tend to broaden the distribution of collec- 
tive linewidths in a 33 x 33 lattice of SRRs. In Fig.f^a), 
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where the lattice spacing is a = 0.5A, the radiative emis- 
sion rates range from the very subradiant 1.2 x 10~ 5 r to 
the superradiant 11T, where T is the decay rate of a single 
meta-atom in isolation. Figure dfb), on the other hand, 
illustrates how the collective effects are diminished when 
the lattice spacing a = 1.4A exceeds a wavelength. The 
distribution of decay rates is considerably narrower with 
the decreased inter-SRR interactions associated with lat- 
tice spacings exceeding a wavelength. Although the ef- 
fects of collective interactions are significantly reduced, 
they do not disappear entirely. The radiative decay rates 
still range from 0.21" to 3I\ 

The dramatically narrowed radiative resonance 
linewidth of some of the collective modes and the 
sensitive dependence of the narrowing on the spatial 
separation of the resonators indicates a strong co- 
operative response of the system to EM fields. For 
very closely-spaced resonators multiple scattering is 
considerably influenced by recurrent scattering events 
in which the field repeatedly scatters from the same 
meta-atoms. In the example studied here, this leads to 
the resonance linewidth narrowing of almost five orders 
of magnitude. Such narrowing could not have been 
described by independent scatterer approach. 

The recurrent scattering that is responsible for the 
dramatic linewidth narrowing can be characterized by 
repeated scattering events between pairs of scatterers, 
triplets of scatterers, etci 20 ' 21 i 23 ~ 25 i 27 In the present work 
we have not analyzed the relative contribution of the dif- 
ferent processes to the distribution of lincwidths. In the 
case of electric dipole scatterers the contribution, for in- 
stance, of repeated exchanges of a photon between pairs 
of dipoles to the distribution of resonance linewidths was 
studied in Ref. ImI . A similar calculation could in princi- 
ple be performed in our system, although the interplay 
between the magnetic and electric dipoles may notably 
complicate the analysis. 

An alternative approach to quantify the contribution 
of different recurrent scattering processes was performed 
in Ref. [25l Numerical simulation results were compared 
with the equations for correlation functions. One, in 
essence, constructs a hierarchy of equations in which the 
nth level describes the recurrent scattering between sub- 
sets of n discrete resonators. Truncating the hierarchy 
after the nth level may therefore be used to quantify the 
contribution of the nth order recurrent scattering. In 
the case of randomly distributed, uncorrclated scatter- 
ers, the role of recurrent scattering between n resonators 
scales with the nth power of density 20 ' 21 i 23 ~ 25 i 27 Correla- 
tions in the positions of the scatterers modify this density 
dependence i 25 1 27 It was found for the both correlated and 
uncorrelated samples 2 ^ that changes in scattering reso- 
nance properties as a function of the density of scatterers 
corresponded to the increased role of recurrent scattering; 
at higher densities the higher order recurrent scattering 
processes become increasingly more important leading to 
the emergence of more strongly subradiant modest 3 - - — 

We now examine the characteristics of some of the col- 




FIG. 3: (color online) An illustration of the most subradiant 
of the collective modes in a 33 x 33 array of SRRs. The height 
of the surface represents the energy of the SRR symmetric 
oscillations |c+| 2 normalized to the peak SRR energy Eq = 
maxf (|c+/| 2 +c_,f| 2 ). The colored patches indicate the phase 
of the electric dipole oscillations for each SRR. The black 
dots indicate the position of each SRR, while their height 
indicates the normalized total energy within each unit cell, so 
that the energy in the magnetic dipole oscillations is given by 
the difference between the black dots' height and the height of 
the surface. This subradiant mode is ferroelectric in nature, 
and has a radiative emission rate 1.2 x 10 _5 r. The lattice 
spacing a = 0.5A, the meta-atom separation within the SRRs 
u = 0.12A, and T E = T M . 

lective modes in a 33 x 33 lattice with an inter-SRR sepa- 
ration of a = 0.5A. As with a single SRR, we can charac- 
terize the state of the system by specifying a complex 
amplitude for both the symmetric (electric) and anti- 
symmetric (magnetic) oscillations. Where the state of 
the system is fully specified by the vector of single meta- 
atom amplitudes (&i, &2, • ■ • , b2N x N y ) T , we represent the 
electric and magnetic oscillations of a single SRR, la- 
belled by £ = 1, ... , N x N y , as c_|_^ and C-^, respectively, 
where 

c±.e = (he-i ± he) • (173) 

As noted earlier, the subwavclcngth proximity of adja- 
cent SRRs permits the creation of extremely subradiant 
collective modes. We illustrate the most subradiant of 
these modes for a lattice spacing of a = 0.5A in Fig. [3] 
The energy of this mode resides almost exclusively in 
symmetric oscillations of the SRRs. However, although 
the meta-atoms in each SRR oscillate symmetrically, the 
electric dipole of each unit-cell resonator element points 
in the opposite direction to that of its nearest neighbor. 
This mode is antifcrroelectric in nature. The phase of 
each electric dipole, indicated by the color of the unit cell, 
forms a checkerboard pattern in the phase profile. This 
mode consists of more strongly excited electric dipole os- 
cillations in the center of the array with smaller contri- 
butions from SRRs on the edges. When this mode is 
excited, the fields emitted from the SRRs tend to remain 
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FIG. 4: (color online) An illustration of the most superra- 
diant of the collective modes in a 33 x 33 array of SRRs. 
The height of the surface represents the energy of the SRR 
antisymmetric oscillations |c_| 2 normalized to the peak SRR 
energy Eo = maxi(\c+,e\ 2 + c_,f| 2 . The black dots indicate 
the position of each SRR, while their height indicates the 
normalized total energy within each unit cell. The total ex- 
citation, |c+,f| 2 , of the electric dipole oscillations is given by 
the difference between the black dots' height and the height 
of the surface. All parameters of the ensemble are as in Fig. 
This superradiant mode has a radiative emission rate of about 
lir and consists largely of magnetic dipole oscillations whose 
phase variation is matched with EM waves propagating in the 
indirections. 



trapped in the ensemble as they repeatedly scatter from 
one meta-atom to another. The scattered fields will leak 
out if this mode very slowly as indicated by the collective 
emission rate of 1.2 x 10" 5 T. 

The most superradiant of the collective modes, shown 
in Fig. [4J by contrast couples very strongly to radiation 
propagating away from the ensemble. This mode is al- 
most entirely magnetic in nature with the SRRs oscillat- 
ing antisymmetrically. These magnetic dipole oscillations 
consist of stripes of constant phase in the y-direction, 
while the phase variation in the a;-dircction is phase 
matched with radiation propagating along ±e x . An EM 
plane wave propagating in the ix-direction whose mag- 
netic field is polarized in the z-direction would have an 
electric field polarized along ±e y . Since the electric 
dipoles in this most superradiant of modes are largely 
unexcited, this mode radiates into an equal superposi- 
tion of EM fields propagating in the positive and neg- 
ative x-directions. The collective excitation coupling to 
these propagating fields results in a spontaneous emission 
rate of HT, more than ten times the single meta-atom 
emission rate. 

In many experimental situations, however, a plane 
wave incident field, with nearly uniform phase and in- 
tensity in the metamaterial plane, drives the ensemble. 
The incident field propagates perpendicular to the plane 
of the metamaterial along the z-direction so that it drives 
the SRRs in phase. It is therefore worthwhile to exam- 
ine modes whose oscillations are phase matched with the 




FIG. 5: (color online) An illustration of the uniform electric 
collective mode in a 33 x 33 array of SRRs. The height of the 
surface represents the excitation energy of the SRR symmetric 
oscillations |c+| 2 normalized to the peak SRR energy Eq — 
maxe (|c+/| 2 + c_^| 2 . The colored patches indicate the phase 
of the electric dipole oscillations for each SRR. The black 
dots indicate the position of each SRR, while their height 
indicates the normalized total energy within each unit cel. 
All parameters of the ensemble are as in Fig. [3] This mode 
consists of the split ring electric dipoles oscillating in phase 
and has a radiative emission rate of approximately the single 
meta-atom emission rate F. 



incident field since they can be addressed directly. The 
two modes of interest are the uniform electric mode, with 
all electric dipoles oscillating in phase, and the uniform 
magnetic mode, where all magnetic dipoles oscillate in 
phase. 

Figure [5] shows the structure of the uniform electric 
mode. As desired, an excitation in this mode has its en- 
ergy almost purely in electric dipole oscillations of the 
split rings. Furthermore, because all electric dipoles os- 
cillate in phase, this mode efficiently couples to EM fields 
propagating out of the plane along ±e z whose electric 
field polarization is along the electric dipoles d = e^. 
Because the fields scattered by this mode propagate out 
of the plane, excitation of the mode by an incident plane 
wave results in reflection of the incident field from the 
metamaterial. In the geometry considered here, the uni- 
form electric mode has a radiative decay rate of 7 e rj T, 
about as strong as the single meta-atom decay rate. 

The second phase matched mode, the uniform mag- 
netic mode, is illustrated in Fig. [6] This uniform mode is 
almost purely magnetic in nature, with all of the meta- 
molecule magnetic moments oscillating in phase, produc- 
ing a sheet of magnetization pointing out of the meta- 
material. In contrast to the uniform electric mode, how- 
ever, this mode cannot strongly couple to fields prop- 
agating out of the plane. In fact, we have found that 
for lattice spacings sufficiently less than a wavelength, 
scattered radiation remains trapped in the ensemble and 
this mode is subradiant. Here, with a lattice spacing of 
a = 0.5 A, the radiative emission rate is suppressed by 
about a factor of 50 below the single meta-atom decay 
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FIG. 6: An illustration of the uniform magnetic collective 
mode in a 33 x 33 array of SRRs. The height of the surface 
represents the energy of the SRR antisymmetric oscillations 
c_ | 2 normalized to the peak SRR energy Eq = max? ( | c+ ^ | 2 + 



the excitation propagates, it begins to broaden so that 
at time t = 500/r, the remaining excitation, containing 
2 x 10~ 4 of the initial energy, has spread through the 
sample. When the lattice spacing is larger, a = 1.5A, the 
excitation spreads more quickly through the sample (at 
time 20/r), indicating weaker EM-mcdiatcd interactions 
between the resonators. In this case only 5 x 10 -8 of the 
initial energy has not been radiated away. 



VII. QUANTIZING THE METAMATERIAL 
DYNAMICS 

In this article, we have developed a general formalism 
to describe collective oscillations in ensembles of meta- 
atoms which comprise a metamaterial. In systems where 
thermal losses are suppressed and can be neglected, how- 
ever, this formalism can easily be quantized. In the quan- 



The black dots indicate the position of each SRR, while tized system, the meta-atom dynamic variables Qj and 



their height indicates the normalized total energy within each 
unit cell. All parameters of the ensemble are as in Fig. [3] This 
mode is composed of all SRRs oscillating antisymmetrically, 
producing magnetic dipoles in phase. The radiative emission 
rate of this mode is 0.02r. 



rate. The form of the magnetic mode does not differ sub- 
stantially from that in Fig. [5] for larger lattice spacings; 
however, inter-resonator spacing affects cooperative in- 
teractions and strongly influence the mode's decay rate^. 
In Ref. it was shown how a subradiant mode analo- 
gous to the uniform magnetic mode we discussed here is 
responsible for the transmission resonance observed in an 
array of ASRsi^ 

The calculated collective modes of the system also de- 
termine the propagation dynamics of localized excita- 
tions. The propagation of excitations are influenced by 
strong interactions between the resonators. Specifically, 
in disordered systems, where the locations of scatterers 
vary randomly, the transition to localization can be char- 
acterized from transport properties^ In the studied sys- 
tem, the positions of the resonators are fixed, so the prop- 
agation dynamics is determined by the particular excita- 
tion. An initial excitation of SRR dipoles will be com- 
prised of some linear combination of collective modes. 
The more radiant components will quickly decay, leav- 
ing behind only the contributions from subradiant modes 
which oscillate at differing frequencies. This behavior 
manifests itself as a decaying propagation and spreading 
of current oscillations through the metamaterial as EM 
fields scatter in the array. The lifetime of the residual 
excitation strongly depends on the presence of recurrent 
scattering and subradiant modes. 

In order to demonstrate the time dynamics of exci- 
tations we have studied the specific example of an ex- 
citation of the left-most strip of magnetic dipoles along 
the y axis in the square array. Such a pattern will lose 
90% of its energy, and propagate a single lattice site in 
a time t = 10 /T for a lattice spacing of a = 0.5 A. As 



their conjugate momenta (f>j, whose Poisson brackets are 
{Qj, 4>j' } = 5j_j>, become quantum mechanical operators 
Qj and 4>j which obey the commutation relations 



Qj,Qj> — 



0j ,(p r 



(174a) 
(174b) 



When quantizing the system, the classical normal vari- 
ables undergo the transformations bj — > s/Kbj and b* — > 

y/hb^. The normal variables thus become harmonic oscil- 
lator creation and annihilation operators which obey the 
commutation relations 



bj,b f 







(175a) 
(175b) 



Similarly, the normal variables for the EM field [Eqs. ([23]) 
and ([23]) ] transform as a q .A -> VMq.A- The EM field 
normal variables then commute with those of the meta- 
atoms and satisfy the commutation relations 



a q ,A,a q ',A' 



^q,A' U q',A' 







3A,A 



<5(q 



(176a) 
(176b) 



The ability to easily quantize this formalism may be use- 
ful in describing the interactions of low loss metamate- 
rials with nonclassical fields. Furthermore, generaliza- 
tions of the formalism to nonlinear metamaterials, e.g., 
involving superconductors, may in and of itself produce 
nonclassical cooperative effects. 



VIII. CONCLUSIONS 

In conclusion, we developed a theoretical formalism to 
describe cooperative interactions of a magnetodiclcctic 
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metamaterial sample with an EM field. We modeled the 
metamatcrial as an ensemble of discrete EM resonators, 
or meta-atoms, that each support a single mode of cur- 
rent oscillation. The meta-atoms could, for example, be 
subwavelcngth circuit elements which support plasmonic 
oscillations. From a Lagrangian describing dynamics of 
the EM field and its interactions with systems of charged 
particles, we derived the conjugate momenta for the EM 
field and meta-atom dynamic variables, as well as Hamil- 
tonian for the metamaterial system. Hamilton's equa- 
tions of motion then describe a coupled dynamics be- 
tween the meta-atoms and the EM field. 

We showed how the EM fields arc emitted from excited 
current oscillations within each meta-atom, and in turn, 
how the EM fields drive the meta-atom dynamics. A 
single meta-atom interacting with its own self-generated 
field behaves as a radiatively damped LC circuit. In an 
ensemble of resonators, the meta-atoms also interact with 
each other. Initially excited by an external field, a meta- 
atom emits EM radiation which then impinges on other 
meta-atoms. The other meta-atoms then re-scatter the 
field. Multiple scattering events mediate an interaction 
between the meta-atoms' current oscillations. The inter- 
actions culminate in a discrete, coupled set of equations 
for the meta-atoms which describe the collective meta- 
material dynamics. The coupled dynamics constituted 
the main results of this article. 

In Sec. |Vl we examined the collective dynamics in 
a regime where the influence of a meta-atom's self- 
generated fields dominates over that of the incident field 
or the fields scattered by all other meta-atoms in the 
metamatcrial. This assumption allowed us to employ the 
rotating wave approximation to simplify the description 
of the dynamics. Appendix ICl on the other hand, gener- 
alized the formalism to provide for a dynamical descrip- 
tion outside the limits of the RWA. 

A metamatcrial possesses as many collective modes as 
there are meta-atoms in the sample, each with its own 
resonance frequency and decay rate. These collective 
modes can behave very differently from oscillations in 
a single, isolated meta-atom. The cooperative interac- 
tions could result in superradiant modes in which energy 
is radiated away more quickly than an ensemble of meta- 
atoms acting independently. Other modes, by contrast, 
are subradiant, for which the mode's radiative emission 
rate is suppressed. As an example, we examined the dy- 
namics of a planar metamaterial formed from a 33 x 33 
square lattice of SRRs. When the resonators are closely 
spaced the collective modes have a broad distribution of 
radiative decay rates. For a lattice spacing of 0.5 A, coop- 
erative interactions suppress the most subradiant mode's 
emission rate by about five orders of magnitude, while 
the most super-radiant mode radiates eleven times faster 
than a single meta-molecule. Finally, we also provided 
an example how the propagation dynamics of excitations 
in a metamaterial array can be analyzed using the collec- 
tive eigenmodes. We found that the lattice spacing, and 
hence the interactions between the resonators, strongly 



influence the rate at which excitations spread over the 
array. In addition to SRRs, the formalism we developed 
could be used to describe interactions between emitters 
with other geometries, e.g. dielectric spheres^. 

The collective dynamics derived from the discrete res- 
onator model can be successfully employed to explain 
experimentally observed phenomena. For example, in 
Ref. H|| we used this model to calculate the resonance 
linewidth narrowing as a function of the system size, as- 
sociated with the experimental observations of the trans- 
mission resonance by Fedotov et al£l. The theoretical 
model provided an excellent agreement with experimen- 
tal findings. This example illustrates how the formal- 
ism developed here lays the ground work allowing one 
to model collective dynamics in large metamatcrial sys- 
tems in which finite-size effects or irregularities may play 
a role. 
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Appendix A: The Lagrangian and the 
Power-Zienau-Woolley transformation 

In this Appendix, we derive the Lagrangian de- 
scribing the dynamics of meta-atoms interacting with 
the EM field given in Eq. (fT3|) . We start from 
the standard Lagrangian for the EM field in the 
Coulomb gauge interacting with arbitrary charge and 
current distributions. Then, using the Power-Zienau- 
Woolley transformation^ - — we express the equivalent 
Lagrangian in terms of polarization and magnetization 
densities. Given the expressions for the polarization and 
magnetization densities in Eq. (|10p . we express the La- 
grangian in terms of effective magnetic fluxes and EMFs 
as in Eq. (fT3| . 

An arbitrary vector field V(r) can be decomposed into 
its longitudinal V|| and transverse Vj_ components 

V = V||+V ± , (Al) 

defined such that V x Vj| =0 and V • = 0. In 
the Coulomb gauge the EM vector potential is set purely 
transverse by requiring that V • A(r) = 0. It follows 
from Maxwell's equations that B(r) is purely transverse 
and that the longitudinal component of the electric field 
E|| is not a true dynamical variable, but is given by an 
algebraic relation by the charge density^ In particular, 
we may write 

By = -VC/, (A2) 
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where 



U{r) = J_ I d * r > 
Aire 



(A3) 



is the scalar potential. The Coulomb energy Vboui is 
given in terms of the meta-atom charge densities in 
Eq. (fT5| . and can be expressed directly in terms of E|| as 



(A4) 



The transverse component of the electric field is given in 
terms of the vector potential A as 



E 



(A5) 



The standard Lagrangian in the Coulomb gauge may 
be written as 



£c — — Vboui + £em + & , 



(A6) 



where 



Ci= j d 3 rj(r,t)-A(r,t) (A7) 

accounts for the interaction between the matter and the 
free EM field, and j = ^ . j 3 is the total current density 
with the contribution from meta-atom j. The meta-atom 
current densities ]j are given in terms of the generalized 
velocities Ij = Qj by Eq. (fT0|) . The vector potential 
A(r, t) provides the continuum of dynamic variables de- 
scribing the evolution of the EM field. The EM field 
dynamics in the absence of charge and current sources is 
governed by the Lagrangian, £ EM [Eq. (fT8|) ]. The charge 
carriers that give rise to the charge and current densities 
have an inertia, and hence the current in a meta-atom, 
resulting from the motion of these carriers, must have 
an associated kinetic energy. This kinetic energy K, is 
given in terms of phenomenological inertial inductances 
in Eq. (fig]). 

The canonical momentum for the fields in the Coulomb 
gauge is given in terms of the time derivative of the vec- 
tor potential and is proportional to the transverse com- 
ponent of the electric field 



n< c )(r) 



dk 



e A(r) = -e E ± (r) 



(A8) 



Similarly, the canonical momentum corresponding to the 
charges Qj is given by 



= ljIj+Xj(t) 



where 



Xj 



J d 3 rA(r,t) ■ [pj(r) + V x Wj -(r)] . 



(A9) 



(A10) 



The factor Xj originates from the interaction Lagrangian 
C\ [Eq. (|A7[) ]: its specific form arises from how the cur- 
rent density jj within each meta-atom j depends on that 



meta-atom's generalized velocity Ij [see Eq. ()10b[) ]. This 
factor represents an averaged projection of the vector po- 
tential onto the current oscillation's mode functions pj 
and Wj. The Hamiltonian in the Coulomb gauge may 
then be derived from the Lagrangian [Eq. (|A6[) ] 



U {C) 



1 

2^ 



[4>i- Xif +?4m + ^Cou1 



(All) 



where the energy of the transverse EM field, or the radi- 
ation field, is responsible for the excitations of the meta- 
atoms 



o/(C) _ e o 



d 3 7 



[Ei(r)+c 2 B 2 (r) 



(A12) 



The quantity Xj originates from the assumption that 
a mode of current oscillation depends on a single dy- 
namic variable with units of charge. The amplitude of the 
charge distribution may change in time, but its spatial 
distribution will not. By contrast, in the more familiar 
scenario where one describes the motion of particles with 
fixed charge qj at a time varying position Tj(t), the con- 
jugate momentum for the position coordinates is given by 
the vector ij + qjA(rj(t)). The scalar quantity \j ar i s ~ 
ing from our model plays the same role as the quantity 
qjA.(vj(t)) appearing in the familiar minimal coupling 
Hamiltonian for moving charged particles. 

Although Eq. (|A11[) is analogous to the standard min- 
imal coupling Hamiltonian description of charged parti- 
cles in an EM field, it does not turn out to be the most 
suitable representation to study the interaction of dis- 
crete scattcrcrs with the EM field. We find it conve- 
nient to express the dynamics in terms of polarization 
and magnetization densities rather than charge and cur- 
rent densities. In this way, when the circuit elements 
are much smaller than a wavelength of EM field with 
which they interact, we may more easily treat the dynam- 
ics in terms of interacting electric and magnetic multi- 
poles. To that end, we employ the Power-Zicnau-Woollcy 
transformation^ For any globally neutral charge distri- 
bution with respective charge and current density p and 
j, there exists a corresponding polarization P and mag- 
netization density M such that 

p(r,t) = -V -P(r,i) (Al3a) 
j(r,t) = P(r,t)+VxM(r,i) (A13b) 

Here, the polarization density is a function of the dy- 
namic variables Qj and the magnetization density is a 
function of their rates of change Ij [Eq. ([9])]. One can 
modify the Lagrangian by adding the total time deriva- 
tive dFJ dt of a function to the original Lagrangian. Here, 
we take 



■ J d 3 rP(r,t) ■ A(r,f) , 



(A14) 



and the equivalent Lagrangian in the length gauge is thus 

dF 



£ = £ 



c 



dt 



(A15) 
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Because F is only a function of the dynamic variables Qj 
and A, the Lagrange equations of motion are invariant 
under this transformation. Explicitly, adding dF/dt to 
the interaction term £7 yields 

£[ = £ I + ^-= j d 3 r (j-p)-A- / d 3 rP-A. (A16) 



From Eq. (|A13bj) . the hrst integral in Eq. (|A16j) can be 

expressed as 



d 3 r (j - P ) • A = / d 3 rA • (V x M) . (A17) 



Integrating this by parts, we obtain the interaction La- 
grangian 



C[ = - j d 3 rB M- d 3 rA • P 



(A18) 



To evaluate the second integral, we recognize that A = 
E + WU, where U(r,t) is the electric scalar potential. 
The last integral appearing in Eq. (|A18|) thus becomes 



d 3 rk ■ P = / d 3 rE • P + / d 3 rV ■ VU. (A19) 



We integrate the last term of Eq. (|A19j) by parts, and 
because U is the Coulomb gauge scalar potential, we ob- 
tain 



J d 3 rP • W = - j d 3 r(V ■ P)U = J d 3 r pU 



2V C 



oul 



(A20) 



Therefore, the Lagrangian in the Power-Zienau-Woolley 
picture can be expressed in terms of the total electric and 
magnetic fields as 



C = K + V c 



oul 



EM 



d 3 r [B • M + E • P] . (A21) 



Although we derived the Lagrangian in Eq. (|A21[) for 
a system composed of ensembles of circuit elements, its 
form is valid for any system of charges where the charge 
density is described by any generalized dynamic vari- 
ables and the current density is a function of their gen- 
eralized velocities. In our system, the total polarization 
P = Ylj P j an( l magnetization M = J2j , with the 
corresponding densities Pj and Mj expressed in terms of 
the dynamic variable Qj and velocity Ij for meta-atom j 
given by Eqs. ([9]). Thus, in an ensemble of meta-atoms, 
the system Lagrangian is given by Eq. (fT3|) . 



Appendix B: Elimination of instantaneous, non-local 
interactions in the Power-Zienau-Woolley picture 

In this Appendix, we provide details of the derivation 
of the Hamiltonian in the length gauge obtained by the 
Power-Zicnau-Woollcy transformation. The derivation is 



analogous to the one discussed in Ref. [50 in determin- 
ing the Power-Zienau-Woolley Hamiltonian for systems 
of charged particles. We begin by examining the portion 
of the Hamiltonian, He [Eq. (|50j) ]. 



He = / d 3 r E ± • D - / d 3 r E • P - V Coul - C 



-EM 



(Bl) 



We will show how the Coulomb potential is absorbed by 
the Power-Zicnau-Woollcy Hamiltonian. We will express 
each term in He in terms of the displacement field D, 
and the polarization density P. The various components 
then combine to yield the Hamiltonian for the free EM 
field, "Hem [Eq. (|3"2"j) ]. the local polarization contact in- 
teraction, and an interaction between the polarization 
density and the displacement field. 

We noted in Sec. IIV Al that an advantage of working 
with the Hamiltonian formalism in the Power-Zienau- 
Woolley picture is that long-range, instantaneous inter- 
actions between meta-atoms do not appear in the Hamil- 
tonian. In any treatment of electrodynamics, the instan- 
taneous non-causal nature of the Coulomb interaction is 
cancelled by other non-causal contributions to dynamics. 
The form of this cancellation, however, is often rather 
subtle. In the Power-Zienau-Woolley Hamiltonian, the 
Coulomb potential is absorbed into a local polarization 
self-energy. Interactions between meta-atoms are then 
mediated entirely by the variables describing the scat- 
tered EM fields. 

In carrying out the simplification, it is useful to note 
the following properties of the longitudinal and trans- 
verse components of any two vector fields Vi and V2. 
The first is that 



d 3 rV 1 ,|(r)-V 2 , ± (r) 



0. 



(B2) 



and as a consequence 
d 3 rV l ■ V 2 



= / d 3 r (VyVj 



V 



2,-L 



(B3) 



We also note, that because the charge density in our en- 
semble of meta-atoms is accounted for entirely by the 
polarization [Eq. (|Al3aj) ]. the displacement field D is 
transverse, i.e., D = D^. We may therefore write the 
transverse and longitudinal electric fields as 

E x = -(D-Px), E||=-1P|,. (B4) 

The Coulomb interaction energy [Eq. (|A4|) ] and the La- 
grangian for the free electromagnetic field [Eq. (|T8j) ] then 
becomes 



^Coul 

£em 



— / 

2e 
1 
2^ 



d 6 r ID 



(B5) 



|B| 2 ) 

i / d 3 r |PJ 2 - — f tfVD • P . (B6) 
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Similarly, the other two integrals appearing in Eq. (|B1[) 
can be expressed as 

J rf 3 rEi • D = J d 3 r (|D| 2 - D ■ p) (B7) 

J rf 3 rE ■ P = i y rf 3 r (d ■ P - |P| 2 ) (B8) 



By the property of Eq. (|B3[) . we may write the portion 
of the Hamiltonian, He as 



H E = n E M + J- I d 3 r IP' 2 



J d 3 r |P| 2 - ^- J d 3 rT> • P, (B9) 



2e 



where Hem, given in Eq. (|32[) . is the Hamiltonian for the 
electromagnetic field. The second term in Eq. (|B9|) has 
absorbed the Coulomb interaction and results only in a 
local meta-atom self- interaction as discussed in Scc. lIV Al 
The final term of Eq. (|B9[) accounts for interaction be- 
tween the distribution of electric dipoles in the polariza- 
tion density and the displacement field. The total Hamil- 
tonian for the system is then given in Eq. (|3ip. 



Appendix C: Collective interactions of strongly 
interacting meta-atoms outside the rotating wave 
approximation 

In Sec. [Vj we saw how the EM field scattered from 
the metamaterial elements produces interactions between 
meta-atom current oscillations in the RWA. For this ap- 
proximation to be strictly valid, the meta-atoms must 
weakly interact with the field, radiatively decaying at 
rates much slower than the oscillator frequencies. In 
many metamaterial systems, however, such assumptions 
can be violated, and the RWA may not be employed. In 
this Appendix, we develop a more general framework for 
the dynamics that allows us to account for very strong 
radiative coupling between the oscillator variables. We 
begin by rcframing the equations of motion for the dy- 
namic variables Qj and their conjugate momenta in terms 
of column vectors of scaled quantities [Eqs. (|135p and 
po) ] 



Q =(Qi,Q2, Qn) t , 



(CI) 
(C2) 



in the frequency domain. For each meta-atom j, the 
scaled charge Qj and its scaled conjugate momentum <j>j 
are slowly varying, with bandwidths comparable to that 
of the incident field's positive frequency component. In 
the time domain, they are related to the physical quan- 
tities Qj and <pj by 



Mt)=yfcC] e-^Q^t) 



(C3) 
(C4) 



In deriving the coupling between the elements, we will 
find that <fi is related to the scaled currents [Eq. ([138|> ] 



I = (Ji,..., In? 



(C5) 



through a dimensionless mutual inductance matrix M, 
and similarly, that the vector of scaled EMFs [Eq. (|116p ] 



£ — (£i, . . . , En) 1 



(C6) 



is related to Q through a matrix resembling a mutual ca- 
pacitance. Since the meta-atoms are separated by signif- 
icant fractions of a wavelength and interactions between 
them are mediated by the radiated field, cj> contains an 
additional contribution from Q. In addition, £ is linearly 
coupled to <fi. This is because oscillating dipoles, whether 
electric or magnetic, produce both electric and magnetic 
fields which drive £ and <f>, respectively. We find that, in 
general, this produces an additional non-trivial coupling 
between resonators. 

We first examine the behavior of the Fourier compo- 
nents of Q and <f> f° r a frequency Q > 0, detuned from 
the central frequency of the incident field by S = Q — SIq. 
Since, by construction, Q, <j), £ and <& are related only 
to the positive frequency components of Qj, <f>j, £j and 
3>j, the scaled variables have no Fourier components for 
6 < — f2o- From the equations of motion for the unsealed 
variables [Eqs. ([38)) ] and the definitions of the scaled vari- 
ables, we arrive at the relations for 5 > — flo 



i{6- 
-i{6 



wt(6) 
£(S) 



(C7a) 
(C7b) 



where u represents a diagonal matrix whose elements 
[<jj]j j = ujj are the resonance frequencies of the indi- 



vidual meta-atoms. Equations (|C7[) represent the cou- 
pling of the meta-atom dynamic variables to the EM 
fields, including the incident field, the fields emitted by 
all other meta-atoms, and the field generated from the 
meta-atom itself. The self-interactions were derived in 
Sec. IV Al [Eqs. ([Mjl and (|95]l]. while we obtained the con- 
tributions from the scattered fields in Sec. IV Bl [Eqs. (|144p 
and (|145p with 6 + £Iq substituted for Oo]- 

As presently written, Eq. (|C7ap . states in terms of 
scaled variables in frequency space that the rate of change 
of the meta-atom charge is equal to its current. Here we 
are interested in how these rates of change are related 
to the states of the meta-atom dynamic variables and 
their conjugate momenta. To express these currents I in 
terms of charges and conjugate momenta, we recognize 
that the conjugate momentum is the sum of the magnetic 
flux and the current multiplied by the kinetic inductance 
lj [Eq. CI])], 



(C8) 



The scaled fluxes, $j [Eq. (|117p ]. contain contributions 
from the meta-atoms' self-generated fields [Eq. (193)) ]. 
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which result in magnetic self-inductances, as well as to 
the magnetic fields generated externally. The kinetic and 
magnetic self-inductances combine to provide the total 
self-inductance [Eq. ([55]) ]. Equation (|139[) thus relates 
the conjugate momentum, <pj , for mcta-atom j to its cur- 
rent Ij and the externally generated flux. The contribu- 
tions to the external flux from the incident magnetic field 
and fields scattered from other meta-atoms in the system 
[Eq. (|145p with fl substituted for Oo] combine to provide 
the external magnetic driving of individual meta-atoms. 
We synthesize these contributions to obtain the relation- 
ship between the column vectors of scaled conjugate mo- 
menta (j>, currents I, charges Q, and fluxes induced by the 
incident field <£> in expressed as 

4>(5) =M(Q.)i(6)-u,- l rl A gT(ri)rl,Q(6)+$ in (5), (C9) 

where we have redefined the diagonal matrices Te and 
Tm containing the meta-atom electric and magnetic 
dipole emission rates so that they reflect the frequency 
dependence of meta-atom scattering rates outside the 
RWA. These matrices have the diagonal matrix elements 

[Te] 3 -j = {^) ( C1 °) 

[Tm],, = (Cll) 
The scaled mutual inductance M, is given by 

M(fl) = (l + iw^TM + w^TigMTi), (C12) 

and the matrices Gm and Gx are given in Eq. (|146j) . The 
diagonal portion of M. has both real and imaginary com- 
ponents: the real part is the self-inductances' contribu- 
tion to this scaled mutual inductance matrix, while the 
imaginary part arises from emission of magnetic dipole 
radiation from the meta-atoms current oscillations. Solv- 
ing Eq. (|C9|) for I, yields 

1(5) =M- 1 (n)(j>(6)-$ in (6) 

+ c- 1 t| i ^(0)T1Q((5)). (C13) 

The current Ij on an clement j is not just related to the 
conjugate momentum <pj, but to the conjugate momenta 
and charges of all other meta-atoms in the system, as 
well as the flux from the incident field. In the absence of 
electric dipole radiation Te = 0, we recover the relation- 
ship between currents and magnetic field fluxes found in 
systems of interacting, radiating, inductive circuits. The 



additional coupling that results from oscillating electric 
dipoles when Te 7^ adds some richness to the dynamics 
of metamaterial systems outside the RWA. 

As with the magnetic fluxes, the EMFs 8 contain con- 
tributions from the self-generated electric fields of the 
meta-atom [Eq. ((94)) ] and from electric fields generated 
by all other meta-atoms in the system [Eq. (|144p with f2 
substituted for f2g]- From the previous results, we can 
express the column vector of EMFs as 

8(5) = - (u - iT E - T10 E (O)T1) Q(5) 

+ T|0 x (fi)T&I(*) +&,(*). (C14) 

The diagonal matrix, u> — iT&, results from the cou- 
pling of each element with its self-generated field where 

T E accounts for decay due to electric dipole radiation. 
j_ 1 

The matrix T E r C/E(^)Y E , where Ge is given in Eq. (|146p . 

provides dipolc-dipolc coupling between electric polariza- 

1 1 

tion densities of distinct meta-atoms, while T|,(7 X (O)T^ 
provides radiated contributions of oscillating magnetic 
dipoles to the EMFs. Since we wish to express the EMFs 
exclusively in terms of the charges and their conjugate 
momenta, we eliminate I by substituting Eq. (|C13p into 
Eq. flUTH) to obtain 

8(5) = - u3- l (Sl)Q(S) + T^x (Q)T M M- X (0)0(5) 

+ 8 in (5)-rig x (n)rl /l M- 1 (n)^ in (5), (ci5) 

where S(f2) is an effective dimensionlcss mutual capaci- 
tance matrix defined such that 

E-\n) = 1 - ilo-'Te - w-^IOb^tI 

- Lo-'riGx (wiM-Htyoj-^gZmri (cie) 

where the matrix expression on the final line arises from 
the expression for current in terms of conjugate momenta 
and charges. This matrix expression is reminiscent of 
a scattering process in which oscillating charges couple 

to oscillating conjugate momenta in other meta-atoms 

1 1 

via T M ^ x (r2)Tg, these conjugate momenta are trans- 
formed into currents by .M (fi), and these currents pro- 
duce electric fields in neighboring meta-atoms through 

Having expressed the currents [Eq. (|C13p ] and EMFs 
[Eq. (|C15|) ] exclusively in terms of charges, conjugate mo- 
menta and driving due to the incident field, we may fi- 
nally write the equations of motion in frequency space for 
the slowly varying charges Q(5) and conjugate momenta 

4,(5). 
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(C17) 

The physical content of Eq. (|C17[) becomes more evident in the limits where the meta-atoms are spaced far enough 
apart that their interactions can be seen as interactions between point sources [i.e., where the approximation in 
Eq. (|147[) holds] and when the magnetic interaction is weak (Tmj <C f2o)- Then, writing Eq. (|C17j) to lowest 
order in Tyi t j/ujj, we obtain dynamic equations where Qj and <pj are driven by fields scattered from the meta-atom 
electric dipoles d^, (fi) = hy y/ujy Cy Qy (S) and magnetic dipoles rn\, (fi) = Ay ujy /Ly ly (5) [sec Eq. (j77j) and 
Eqs. P5]l . ip5j) . ipgj) and (fWl) ] 

-iflQ 3 (S) H 



Hok 



E(°< 

TT3 



&(<5) 



cGx^-ry.^.d^Cn) 



(C18a) 



-iQ,(f>j(5) 



I — 1 



Qi(<y) = 6 



in,J 



J2(G(r j -ry,n).d ( y + \Q) + -G x (r 



iy,fi). m «(Q) 



(C18b) 



where otj = (2Aj /2>)yjujj/ Lj /3 and Xj = ^hj^/ojjCj/i. The incident field drives the meta-atom electric and magnetic 
dipoles, producing the terms — ujy^ m j and £i n ,j- The effects of the meta-atoms' self-generated fields are included 
on the left hand side of Eqs. (|C18|) . The magnetic fields produced by all other meta-atomic dipoles in the system 
j' j drive the dynamics of Qj, while the electric fields scattered from these meta-atoms drive the dynamics of the 
conjugate momentum <fij. 



In principle, one could solve Eq. ()C17[) in the narrow 
bandwidth approximation in which the driving field en- 
velopes Ej n and B; n vary on time scales much larger than 
l/fio. O ne would accomplish this by substituting Oo for 

in the coupling matrices Ge/ m/x a n d T E / M , then in- 
verse Fourier transforming Eq. ([Crf) . This proced ure, 
however, may not be particularly illuminating. We find 
it useful to explore the dynamics in terms of the oscilla- 
tor normal variables b. But first, we will revisit the basic 
characteristics of these normal variables to understand 
how they behave outside the RWA. 

1. The normal meta-atom variables revisited 

In Sec. |Vl we have assumed the meta-atom normal 
variables b = (&i, 62, . . . , 6at) t are slowly varying, and 
are proportional to the slowly varying envelopes of the 
charges Q and conjugate momenta <f>. This was a good 
approximation when we assumed each meta-atom cou- 
pled to its self-generated fields much more strongly than 
it couples to the external fields, i.e., when Tsj, ^mj "C 
Wj, f2o- in those limits, the external field interactions act 
as a perturbation that slowly alters the dominant behav- 
ior of the meta-atoms oscillating as effective LC circuits. 
When the coupling to the external field is stronger, how- 



ever, we will see that this is no longer the case. 

While the normal variables are slowly varying in the 
RWA, in general they contain fast oscillating components 
even when Q and <f> have narrow bandwidths. To see this, 
we rewrite the normal variables [defined in Eq. f| 109[) ] in 
terms of the slowly varying dynamic variables. Recall 
that in terms of the slowly varying quantities, the original 
charges and conjugate momenta are expressed as 

Q 1^L = e- inot Qy{t) + e inot Q*(t) (C19a) 
J^£L = e-Mot^.^ + e «io*0»( t ). (C19b) 

We therefore express the column vector of normal vari- 
ables as 

b(i) - bW(i) + e 2tn «[b^{t)}* , (C20) 

where we have defined the slowly varying normal variable 
components 

bW(i)=^(Q(i)±#)). (C21) 

The normal variables b therefore have a bimodal spec- 
trum with a slowly varying component peaked at zero 
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frequency and fast oscillating component peaked at fre- menta from Eq. (|C17[) . In terms of the normal variables 
quency — 2f2 . Note that despite the a ppar ent similar- the vectors Q and 4> ar e given by 
ity between the definition of bj [Eq. ()109p ] and b^ 



b(+) 



In the 



[Eq. (|C2ip ]. b< _ ' is not equal to 

RWA, for example, we ignore the rapidly rotating term 
[cxp(2ifi )[b(~)(£)]* in Eq. IC20| and make the approxi- 
mation b^~' w 0. Physically, the RWA implies that the 

normalized conjugate momentum (t) has the same am- / q\ 1 / 1 1 \ / fj( + ' 

plitude as Qj(t) but its oscillation lags by a definite phase I ^ J ~~ \ —i i J I b(~) 

7r/2 in accordance with Eq. (|137[) . However, outside the 
RWA, the contributions of b 1 - - ' cannot be neglected. 



(C22) 



2. Normal variable dynamics outside the RWA 

We can obtain the normal variable dynamics from The vectors of normal variables b^) therefore evolve ac- 
those for the slowly varying charges and conjugate mo- cording to 



~ [ U-)(s) ) - { c<-.+)(n) c(-.-)(n) ) { £<->(*) ) + [ 4-) (5) ) ' (C23) 

where C (si ^ (si,s 2 = ±) are N x N matrices that provide a linear coupling between the column vectors b^ 2 ^ and 
t>( Sl '. These coupling matrices are given by 

Ms, s,1 ■ /o \r SlT E + S2T M A^ _1 

C (Sl ' S2j =t(Q - siu) S SuS2 



2 



(C24) 



The top line of Eq. (|C24p contains the diagonal elements of the matrices C' ai '* 2 ^ that arise from the interaction of 
the meta-atoms with their self-generated fields. For example the diagonal elements of C^ +,+ - ) contain the detunings 
A [Eq. (|152p ] and the total radiative decay rate Ye + Tm- The interaction matrices C^ S1,S2 ' contain the effects of all 
scattering processes, including those resulting from scattered electric fields emitted from electric and magnetic dipolcs 
and those resulting from magnetic fields emitted from magnetic and electric dipolcs. Interaction with the incident 
field produces the driving represented by the column vectors 



f£\6) = -j= [±iS - (wizTgSxT^J M- 1 ^ (C25) 

The EM mediated interactions simplify greatly if we assume the self-inductance of a meta-atom is much greater 
than the mutual inductance between any two meta-atoms. A necessary condition for this is that Tyi,j *C uij for all 
j since M = 1 + OiTuUJ' 1 ). In this limit, we neglect all contributions of order Tyij/uij to the mutual inductance, 
allowing us to make the substitution Ai^ 1 ss 1. This yields 

C<i'«> « i (flo - biw) S SUS2 ' SiTe + S2Tm + ^n^n + ^Tigurj + TjggTj+ jif2 T|g^ 



f 



Under the additional assumption that all meta-atom resonance frequencies lie in a narrow bandwidth around flo, the 
matrix providing the dynamic coupling between the vario 
normal variables in the RWA [Eq. ([T5T])]. i.e. C(+>+) w C. 



matrix providing the dynamic coupling between the various bj and bj, is identical to the coupling matrix between 
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3. Temporal dynamics and collective modes outside the RWA 

When the incident field possesses a narrow bandwidth around £1q and varies much more slowly than the time it 
takes for light to propagate across the metamaterial sample, we can obtain a simple expression for the collective 
temporal evolution of the metamaterial. With this slowly varying incident field, we can approximate the dynamics 
by replacing the frequencies £1 = Qq + 5 appearing in the interaction matrices with D,q. We then inverse Fourier 
transform Eq. (|C23j) to obtain 



d 



b<+>(t) 
b(-)(t) 



C<+'+>(fio) C(+--)(n ) 
C(-'+)(O ) C(-'-)(O ) 



b(+)(t) 
b(~)(t) 



£->(*) 



(C27) 



These equations describe the collective response of a metamaterial to a narrow bandwidth incident field where the 
inter-resonator interactions and emission rates can be arbitrarily large. Unlike the simplified collective dynamics 
derived in Sec. IV Bl Eq. (|C27|) is not subject to the constraints of the RWA. 



4. Recovering the dynamics of the RWA 

The dynamics in the RWA that we explored earlier in 
Sec. fVl amounted to neglecting the fast oscillating com- 
ponents of the normal variables b. Here, this equates 
to assuming b^ - ) = 0, and therefore b = t/ + ). We ar- 
gued earlier that this approximation is valid in the lim- 
its of weak interaction - i.e., Tej, Tmj <C ilo ~ and 
in which all single meta-atom resonance frequencies lie 
within a narrow bandwidth about the driving frequency 
- i.e., |Aj| -c f2o. Indeed, when the interactions are suffi- 
ciently weak, the diagonal elements of C' _,_ ) [Eq. (|C23[) ] 
[i(O + lj) + (Te + Tm)/2] dominate over every other ele- 
ment of the coupling. As a result, in the response of the 



metamaterial to the incident field, the elements of b^ - 1 
would be negligible in comparison to t>M. Note, how- 
ever, that although Te/m <C is a necessary condition 
for the validity of the RWA, it is not sufficient in and of it- 
self. This is because the interaction between elements can 
still become very strong if the separation between them 
is much less than a wavelength. Here, however, we will 
assume the inter-element separation is sufficiently large 
that the limits on T E / M are sufficient. In the RWA, we 
may therefore expand C( + > + ) to lowest order in Te and 
Tm and make the approximation (ft/cuj) 3 w 1. Upon 
doing this, we recover precisely the dynamics given in 
Eq. (fl"49|) of subsection |VB] 
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